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SECTION  1 
INTRODUCTION 

Linear  fracture  mechanics  has  gained  a substantial  ac- 


ceptance in  industries  and  has  become  one  of  the  most  impor- 
tant design  considerations.  It  is  now  well  recognized  that 
many  engineering  structures  such  as  airplanes,  rockets,  pres- 
sure vessels T bridges,  etc.  may  contain  preexisting  flaws  and 
that  crack  propagation  emanating  from  such  flaws  occupies  a 
major  fraction  of  useful  fatigue  life  of  structures.  These 
structures  are  being  designed  for  fracture  prevention  and  for 
limited  service  life.  The  basic  parameters  in  linear  fracture 
mechanics  are  the  elastic  stress  intensity  factors  Kj,  and 

KII;[  for  the  opening  mode,  the  in-plane  shear  mode  and  the 


anti-plane  shear  mode,  respectively.  Numerous  methods  are  now 
available  for  determining  these  factors. 

One  of  the  most  powerful  methods  which  can  be  used  to  de- 
termine the  stress  intensity  factors  in  structures  of  compli- 
cated geometry  is  the  finite  element  method.  A large  number 
of  papers,  most  of  which  have  appeared  since  1970,  discuss  the 
application  of  the  finite  element  methods  to  fracture  mechan- 
ics. The  survey  on  the  finite  element  analyses  of  two-dimen-  * 

sional  fracture  problems  has  been  made  by  several  authors  [1- 
7]  . 

•fj  When  applying  the  finite  element  method  to  analyze  frac- 

ture problems  which  contain  regions  of  stress  singularities. 
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it  is  necessary  to  increase  the  detail  of  the  finite  element 
model  in  the  vicinity  of  the  singularity  to  obtain  reasonably 
accurate  results.  However,  most  attempts  to  model  a crack 
front  by  a profusion  of  tiny  conventional  elements  have  been 
unsuccessful  due  to  the  limited  storage  capacity  of  the  pres- 
ent-generation computer  and  to  the  resulting  malconditioned 
global  stiffness  matrix.  In  reference  8,  it  has  been  shown 
that  the  convergence  rate  for  the  finite  element  solution  in 
the  presence  of  singularities  is  controlled  by  the  nature  of 
singularities  and  that  the  conventional  high  accuracy  element 
obtained  by  using  polynomials  of  high  order  as  interpolation 
functions  cannot  improve  the  rate  of  convergence  and  they  are 
either  inadequate  or  inefficient. 

During  the  last  five  years,  to  overcome  such  difficulty, 
there  have  been  many  investigators  who  developed  special  crack 
elements  which  incorporate  the  singular  behavior.  These  spe- 
cial elements  are  used  to  model  the  region  near  the  crack 
front,  while  the  remaining  region,  in  general,  is  modelled 
with  conventional  elements.  Owing  to  these  special  crack  ele- 
ments, the  linear  elastic,  two-dimensional  fracture  problems 
seem  to  be  well  in  hand  for  application  of  the  finite  element 

I 

methods . ; 

Although  there  has  been  a notable  increase  in  interest  in 
the  application  of  the  linear  fracture  mechanics  to  realistic 
engineering  problems  which  are  inherently  three-dimensional, 
only  a few  studies  in  finite  element  analysis  of  such  problems 
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have  been  published  to  date  [56,58,59,63-70,72,75,76].  Par- 
ticularly, for  the  problem  of  surface  flaw  which  is  the  most 
common  type  of  realistic  crack  no  analytical  solution  has  been 
reported  to  date.  Numerical-analytical  method  such  as  an  al- 
ternating technique  has  been  limited  to  problems  of  simple  ge- 
ometries [39,40,42,43,45].  Huge  amount  of  computation  time 
and  storage  requirement  of  matrices  of  large  bandwidths  far- 
ther restrict  the  use  of  the  conventional  finite  element  anal- 
ysis. 

In  three-dimensional  solid,  the  stress  in  the  plane  nor- 
mal to  and  near  the  crack  front  has  the  same  angular  distribution 
as  that  of  a two-dimensional  crack  problem  under  the  action  of 
in-plane  stretching  and  shear  and  out-of-plane  shear  [38].  It 
has  an  inverse  square  root  r singularity.  In  the  last  few 
years  several  authors  have  developed  three-dimensional  crack 
front  elements,  taking  such  kind  of  singularity  at  the  crack 
front  into  consideration  [58,59,66,67,72,75].  However,  effort 
is  still  needed  to  seek  for  efficient  methods  to  obtain  reli- 
able stress  intensity  factor  solutions. 

Another  important  crack  problem  is  to  determine  the  plate 
bending  and  shearing  stress  intensity  factors,  and  Kg,  for 
a sharp  notch  in  a plate  subjected  to  out  of  plane  bending. 
There  have  been  several  analytical  studies  dealing  with  the 
bending  of  a cracked  plate.  Most  recently,  Ishida  analyzed 
the  interference  effect  of  cracks  located  randomly  in  an  in- 
finite plate  subjected  to  out  of  plane  bending  based  on  the 
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Poisson-Kirchhoff  thin  plate  theory. 


* 

? 

The  assumed  stress  hybrid  finite  element  method  pioneered 
by  Pian  [85],  has  been  successfully  applied  to  problems  in 
which  conventional  finite  element  model  finds  difficulties. 

f 

Notable  examples  are  bending  of  thin  plates  and/or  thick 
plates,  problems  with  prescribed  boundary  conditions,  shell 
problems  and  singularity  problems. 

I • There  exist  two  assumed  stress  hybrid  finite  element  mod- 

els for  the  development  of  special  crack  elements  for  linear 
fracture  mechanics:  (1)  Scheme  I is  based  on  the  complete  ex- 
pansion of  the  exact  near-field  solution  each  term  of  which 

\ • 

satisfies  both  equilibrium  and  compatibility  conditions 
- inside  the  element.  (2)  Scheme  II  is  based  on  assumed  equi- 

librating  stress  field  including  the  singular  term  and  on  in- 
dependently assumed  boundary  displacements.  The  singular 
stress  term  also  satisfies  the  compatibility  condition. 

In  the  first  part  of  this  thesis,  special  crack  front 
elements  are  developed  for  direct  evaluation  of  stress  inten- 
sity factors  Kj,  Kjj  and  KII3;  of  arbitrarily  shaped  three-di- 
mensional cracks  based  on  the  second  scheme  of  the  hybrid 
stress  finite  element  model  described  above.  The  elements 
, have  been  utilized  to  analyze  several  fracture  test  specimens 

^ of  typical  geometries,  i.e.,  a center  crack  tension  specimen, 

a single  edge  crack  tension  specimen,  a double  edge  crack  ten- 
M 

'l  sion  specimen  and  a compact  tension  specimen.  Other  numerical 

examples  include  an  axisymmetrically  located  penny-shaped 


crack  in  a cylindrical  rod,  an  embedded  penny-shaped  crack  in 
a cube,  a semi-circular  surface  crack  in  a rectangular  prism 
and  a quarter-circular  corner  crack  in  a rectangular  prism 
under  uniform  tension  loading. 

In  the  second  part  of  this  thesis,  a superelement  is  de- 
veloped for  the  analysis  of  plate  bending  and  plate  shearing 
stress  intensity  factors.  Kg  and  Kg,  of  thin  plates  with  a 
through-the- thickness  crack  subjected  to  out  of  plane  bending. 
The  particular  approach  is  based  on  the  first  scheme  of  the 
hybrid  finite  element  model,  Poisson-Kirchhoff  thin  plate 
theory  and  the  complex  variable  technique.  Numerical  example 
includes  the  problem  of  pure  cylindrical  bending  of  a thin 
plate  with  finite  width  W and  a centrally  located  through-the- 
thickness  crack  of  length  2a  for  a number  of  cases  with  dif- 
ferent ratios  of  2a/W.  This  method  of  approach  is  extended  to 
bending  analysis  of  a bi-material  plate  with  through-the- 
thickness  cracks  located  along  and/or  normal  to  the  interface, 
an  isotropic  plate  with  a wedge-shaped  notch  and  an  anisotrop- 
ic plate  with  through- the- thickness  cracks. 
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SECTION  2 


ON  LINEAR  FRACTURE  MECHANICS 
2.1  Introduction 

By  considering  the  configuration  of  an  isolated  slit 
crack  which  minimizes  the  total  free  energy  of  the  system, 
Griffith  first  formulated  a basic  criterion  for  crack  exten- 
sion in  terms  of  the  fundamental  energy-balance  concept  of 
classical  mechanics  and  thermodynamics  [29].  The  principles 
laid  down  in  his  pioneering  work  and  implications  drawn  from 
them  have  played  a vital  part  in  the  development  of  present 
day  fracture  mechanics  theory.  The  Griffith's  criterion  can 
be  written  as 


acr  = ( 2Ey/Tra)  ** 

for 

plane 

stress 

condition 

and 

acr  = { 2Ey/ira  (1  - v*)}*5 

for 

plane 

strain 

condition 

where  a is  the  critical  applied  stress,  E is  the  Young's 

C* 

modulus,  y is  the  free  surface  energy  per  unit  area,  a is  the 

half  length  of  the  crack  and  v is  the  Poisson's  ratio.  If  the 

applied  stress  exceeds  the  critical  value  of  o , a brittle 

cr 

crack  is  free  to  propagate  spontaneously  without  limit. 

Irwin  gave  the  Griffith's  concept  a more  general  theoret- 
ical framework  [31].  In  his  theory  the  concept  of  the  elastic 
stress  intensity  factor  K which  is  the  magnitude  of  the  local 
stress  field  is  introduced  in  place  of  the  Griffith's  free 
surface  energy  y . It  can  be  shown  that  both  the  energy  ap- 


proach  and  the  stress  intensity  factor  approach  are  equivalent 
and  the  stress  intensity  factor  is  related  to  the  free  surface 

energy  through  the  strain  energy  release  rate  G = (9U/9a)  as 

j*  2 k 

K = (EG)  for  plane  stress  problems  and  as  K = {EG/(1  - v ) } 

for  plane  strain  problems.  The  Griff ith-Irwin  fracture  crite- 
rion states  that  the  fast  fracture  occurs  when  K or  G exceeds 

a critical  value  K or  G . K and  G are  material  constants 

c c c c 

and  are  refered  to  as  the  fracture  toughness  of  the  material. 
They  have  the  dimension  of  stress* (length)  and  stress • length 
respectively. 

2.2  Modes  of  Fracture 

Irwin  also  proposed  three  basic  modes  of  crack  surface 
displacement  as  in  Figure  1.  Mode  I (opening  mode)  corre- 
sponds to  normal  separation  of  the  crack  surfaces  under  the 
action  of  tensile  stresses;  mode  II  (sliding  mode)  corresponds 
to  mutual  shearing  of  the  crack  surfaces  in  a direction  normal 
to  the  crack  front;  mode  III  (tearing  mode)  corresponds  to  mu- 
tual shearing  parallel  to  the  crack  front.  Any  crack  problems 
can  be  considered  to  be  a combination  of  these  modes  for  var- 
ious loading  situations.  However,  since  there  is  always  a 
tendency  for  a brittle  crack  to  propagate  in  the  direction 
which  minimizes  the  shear  loading,  the  first  mode  is  by  far 
the  most  pertinent  to  crack  growth  in  brittle  materials.  In 
fact,  to  the  author's  knowledge,  the  experimental  value  of  Gc 

or  K for  mode  II  and  mode  III  crack  extensions  has  not  been 
c 

reported  to  date. 
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2.3  Asymptotic  Field  in  the  Vicinity  of  the  Crack  Front 

Sneddon  first  obtained  the  local  stress  field  in  the  vi- 
cinity of  the  crack  front  for  a penny- shaped  crack  [32] . Rec- 
ognizing the  general  applicability  of  this  stress  field, 
Williams  extended. the  inverse  square  root  r character  of  the 
local  stresses  of  Sneddon's  solution  to  the  most  general  case 
of  two-dimensional  isotropic  crack  problems  [30].  Later, 

Irwin  reformulated  the  local  stress  and  displacement  field 
based  on  his  characterization  of  three  basic  fracture  modes 
[31].  In  their  simplified  crack  model,  it  is  assumed  that  the 
crack  f?ont  is  perfectly  sharp  and  that  the  crack  surface 
remains  stress  free  at  all  stages  of  loading.  The  stress  and 
displacement  fields  thus  obtained  take  a simple  analytic  form 
and  are  given  below  for  each  of  the  three  modes  of  fracture 
under  plane-strain  condition: 

Mode  I: 

ox  = {Kj/ (2irr)  ^Icosl  (1  - sinlsin3^) 
ay  = v(ax  + az> 

a = {K  /(2TTr)J5}cos|(l  + sin^sin3^9 ) 

Z I 2 2.  I 

Tzx  = {Kj/ (27Tr)  Js}sin|cos|cos3^9 

T - T - 0 
xy  yz 

u = (Kj/G)  (r/2ir ) ^cos^  (1  - 2v  + sin2|) 
v = 0 

w = (Kj/G)  (r/2ir)  ^sinl  (2  - 2v  - cos2|) 
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Mode  II: 


ox  = -{KII/ (2-nrr)  ij}sin|  (2  + cos|cos3^0) 
ay  = v(ax  + V 

oz  = {KI;I/ (2irr)  35}sin|cos|cos329 

tzx  = {KXj/  (2Trr)  ^}cos|  (1  - sin|sin3^9) 

t — t =0 
xy  yz 

u ~ (r/27r)  ^sinl  (2  - 2v  + cos2|) 

v * 0 

w = (K^/G)  (r/2Tr)  ^cosl  (-1  + 2v  + sin2|) 

Mode  Ills 

a 2 a - a = t =0 
x y z zx 

Txy  * ~{KIII/(2^r),5}sin| 

Tyz  = {KIII/(27rr)35}cosl 
u = w = 0 

v = (Ki;ei/G)  (2r/ir)  ^sinl 

where  the  crack  lies  along  the  x-axis' in  the  y-plane  and  the 
polar  coordinates  (r,0)  are  measured  from  a crack  front 
[Figure  2],  G is  the  elastic  shear  modulus,  v is  the  Poisson's 
ratio  and  the  Kj,  and  KII];  are  stress  intensity  factors 

corresponding  to  mode  I,  mode  II  and  mode  III  type  of  fracture 
respectively.  In  the  asymptotic  stress  and  displacement  field 
described  above,  the  stress  intensity  factors  which  specify 
the  magnitude  of  the  local  field  depend  only  on  the  applied 
load  and  the  geometry  of  the  cracked  body.  The  remaining 
terms,  which  are  divided  further  into  a radial  component  and 
an  angular  component,  depend  only  on  spatial  coordinates  about 


1 


sional  crack  under  plane  strain  condition,  i.e., 

ox  - (l//2-rrr)  {Kjcosl  (1  - sin|sin3^)  “ K^sin^^  + coslcos3^)} 

ay  = (2v//2irr)  {Kjcosl  - K^sin^ 

az  = (l//27rr)  {Kjcosl  (1  + sin|sin3^9)  - Kjjsinlcoslcos3^ } 

t = (1/ /2irr)  {K  sintcosfcos3^  + K__cos|(l  - sin^sin\9)} 

XZ  1 2 2 2 II  2 2 2 

Tyz  = (l//2TTr)KIIlCos| 

Txy  = -(l//2ur)KII3;sin| 

u = (1/8G)  (2r/iT)  [ (5  - 8v)cos|  - cos3|]  + K^lO  - 8v)  sin 

+ sin3^9  ] } 

» 

v = (1/G)  (2r/7r)  ^KjjjSinl 

w =*  (1/8G)  (2r/ir)  [ (7  - 8v) sin|  - sin3^9  ] - K^tO  “ 8v)cos 

+ cos3^9  ] } 
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Fracture  Analysis 

In  this  section,  a brief  review  is  given  of  numerical 
techniques  other  than  finite  element  method  that  are  readily 
applied  to  three-dimensional  linear  fracture  problems.  The 
finite  difference  method  is  omitted  because  its  application  to 
three-dimensional  crack  problems  has  been  limited  to  elastic- 
perfectly  plastic  analysis  [41]. 

(1)  Schwartz-Newmann  Alternating  Technique 

The  Schwartz-Newmann  alternating  technique  is  an  itera- 
tive method  [42] . In  using  this  method  stresses  are  first  ap- 
plied to  the  crack  surface  and  using  the  known  analytical  so- 
lution for  an  embedded  crack  stresses  can  be  computed  in  a 
semi-infinite  solid  at  the  locations  of  the  surfaces  of  the 
solid  in  which  the  crack  is  assumed  to  reside.  Then,  the  op- 
posite of  these  stresses  is  applied  to  the  surface  of  the  body 
in  question.  The  resulting  stresses  at  the  location  of  the 
crack  surface  are  computed  by  using  the  known  analytical  solu- 
tion of  a semi-infinite  solid  subjected  to  variable  surface 
traction  over  its  flat  surface.  The  crack  surface  stresses 
are  again  removed  as  before  and  the  process  is  repeated  until 
further  contributions  are  negligible.  The  results  of  all  it- 


eration are  superimposed  to  obtain  stress  intensity  factors. 


This  method  requires  an  analytical  solution  for  an  infinite 
solid  with  the  crack  subjected  to  variable  pressure.  Since 
such  solution  is  not  available  for  cracks  of  complex  geome- 
tries, the  alternating  technique  can  be  applied  only  to  cracks 
of  simple  geometry. 

(2)  Boundary  Integral  Equation  Method 

Boundary  integral  equation  method  has  been  used  exten- 
sively to  solve  practical  three-dimensional  crack  problems  by 
Cruse  et  al.  [46,47,50].  This  method  also  uses  the  known  sin- 
gular solution  to  the  Navier-Cauchy  elasticity  equations  and 
this  gives  rise  to  a vector  identity  similar  to  that  of 
Green's  formula  for  Laplace's  equation.  Taking  the  reference 
area  to  lie  on  the  boundary  of  the  body  in  question,  a set  of 
simultaneous  integral  equations  which  relates  boundary  dis- 
placements to  corresponding  boundary  tractions  is  generated. 
Since  either  of  these  boundary  quantities  determines  the  other 
the  unknown  displacements  or  tractions  of  each  boundary  refer- 
ence area  are  obtained  by  solving  the  resulting  simultaneous 
equations.  The  boundary  quantities  are  assumed  to  be  constant 
or  linear  over  each  reference  area. 

2.5  Finite  Element  Methods  for  Three-Dimensional  crack 

Problems  !. 

2.5.1  Introduction 

Finite  element  methods  have  been  used  in  linear  fracture 
mechanics  quite  extensively.  The  survey  on  the  finite  element 
analyses  for  two-dimensional  fracture  problems  has  been  made 


by  several  authors  [1-7] . 

The  application  of  finite  element  methods  to  three-dimen- 
sional linear  fracture  problems  can  be  classified  into  the 
following  three  categories: 

(1)  Conventional  finite  element 

(2)  Special  singular  element 

(3)  Superposition  of  analytical  and  finite  element  solutions 

2.5.2  Conventional  Finite  Element  Model 
The  assumed  displacement  approach,  in  general,  cannot 
provide  the  stress  intensity  factor  solution  directly. 

However,  various  quantities,  such  as  displacements,  stresses, 
strain  energy  release  rate  and  local  strain  energy  are  relat- 
ed to  the  stress  intensity  factors.  Therefore,  once  one  of 
these  quantities  is  obtained  by  the  finite  element  analysis, 
the  stress  intensity  factor  can  be  calculated.  Various  proce- 
dures have  been  developed  to  translate  the  finite  element  so- 
lutions to  the  stress  intensity  factors.  Basically,  they  are 
divided  into  the  following  three  categories. 

(1)  Based  on  a crack  opening  displacement  (C.O.D.) 

The  mode  I stress  intensity  factor  can  be  represented  in 
terms  of  local  coordinates  (x,y,z)  and  (r,0)  as  in  Equation 
(2.2).  Thus  knowing  w at  some  point  on  the  crack  surface,  the 
stress  intensity  factor  K ^ is  obtained  by  substituting  w into 
Equation  (2.2).  More  accurate  K^.  solution  can  be  obtained  by 
extrapolating  the  solutions  of  several  reference  points  to  the 


crack  front. 


(2)  Based  on  a stress  value 

Normal  stress  oz  is  related  to  the  mode  I stress  intensity 

factor  Kt  by  Equation  (2.1).  Thus,  if  the  value  of  a at  some 

point  near  the  crack  front  is  obtained  £rom  the  finite  element 

/ 

analysis,  Kj  is  calculated  readily.  In  this  case,  extrapola- 
tion technique  can  also  be  used  to  obtain  accurate  Kj  solution. 

(3)  Based  on  a strain  energy  release  rate 

The  stress  intensity  factor  Kj  is  related  to  the  strain 
energy  release  rate  by  the  following  equation. 

G = K2 (l-v2)/E 

First  the  total  strain  energy  for  a given  crack  geometry  is 
computed  and  then  the  total  strain  energy  for  a small  crack 
extension  is  recomputed.  The  resulting  total  strain  energy 
difference  AU  divided  by  the  difference  in  surface  areas  of 
the  two  crack  lengths  AA  yields  the  strain  energy  release 
rate  G as 

r = Slim  AU 
~ Aa-*o  AA 

However,  this  procedure  cannot  be  used  effectively  in  three- 
dimensional  analysis  br cause  of  the  variability  of  KT  along 

» i. 

the  crack  front. 

Miyamoto  and  Miyoshi  [55-57]  analyzed  a semi-elliptical 
surface  crack  in  a plate  of  finite  thickness  subjected  to 
tension  by  conventional  finite  element  method.  They  used  six 
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node-pentahedron  elements  built-up  from  three  four- node 
tetrahedron  elements.  The  crack  opening  displacement  neces- 
sary for  the  determination  of  stress  intensity  factor  Kj  was 
obtained  through  a two  step  computation  involving  coarse  mesh 
(2332  d.o.f.)  analysis  for  the  entire  plate  and  fine  mesh 
(1722  d.o.f.)  analysis  for  the  small  region  near  the  crack 
front.  The  reported  CPU  time  is  1814  seconds  -for  the  step  1 
and  1457  seconds  for  the  step  2 analysis. 

Bergan  and  Aamodt  [62]  employed  a conventional  twenty 
node  isoparametric  hexahedron  elements  to  study  the  crack 
propagation  of  a semi-elliptical  surface  crack  in  a thick 
plate.  For  increasing  the  computational  efficiency,  they 
utilized  certain  techniques  such  as  multilevel  subdivision 
of  the  structure  (substructure)  and  condensation  of  the 
unnecessary  degrees  of  freedom.  The  stress  intensity  factor 
was  calculated  based  on  the  strain  energy  release  rate  G. 

Three  levels  substructures  contained  3360  nodal  degrees  of 
freedom.  For  obtaining  the  condensed  stiffness  matrix  asso- 
ciated with  360  nodal  degrees  of  freedom,  a computing  time  of 
one  and  a half  hours  were  required,  while  only  1.6  seconds 
of  CPU  time  were  required  to  solve  the  condensed  matrix  equa- 
tion. 

Miyata,  Shida  and  Kusumoto  [64]  developed  a simple  method 
for  estimating  stress  intensity  factors  by  using  only  four- 
node  constant  strain  tetrahedron  elements  near  the  crack 


front.  They  obtained  a simple  formula  to  relate  the  stress 
intensity  factors  to  the  average  stresses  in  a given  area 
around  the  crack  front.  Initially  this  basic  approach  was 
developed  for  two-dimensional  crack  analysis  by  using  constant 
strain  triangular  elements  at  the  crack  tip  and  later  extended 
to  three-dimensional  analysis.  However,  even  the  mesh  sub- 
division for  the  two-dimensional  analysis  of  the  center  crack 
specimen  contained  as  many  as  several  thousand  degrees  of 
freedom. 

Parks  presented  a technique  to  obtain  stress  intensity 
factors  by  a single  conventional  finite  element  analysis 
based  on  the  energy  release  rate  [61] . In  his  method,  the 
energy  release  rate  is  expressed  as 


8 tt  .2 

e = k2 
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If  the  displacement  q and  the  rate  of  changing  stiffness 
matrix  K are  determined,  the  stress  intensity  factor  is 
readily  obtained  from  above  equation.  This  procedure  was 
extended  by  Parks  [61]  to  general  three-dimensional  crack 
configurations  and  was  applied  to  a three-dimensional  formula- 
tion of  an  axisymmetric  penny  shape  crack  problem. 

2.5.3  Special  Singularity  Element 

Although  the  complexities  of  formulating  three-dimen- 
sional singularity  element  are  multiplied  manyfold  over  those 
of  analogous  two-dimensional  element,  several  authors 
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developed  three-dimensional  crack  front  elements,  taking  the 
singularity  at  the  crack  front  into  consideration. 

Based  on  assumed  displacement  finite  element  model, 

Tracey  [58,59]  developed  six- node  pentahedron  elements  for 
the  three-dimensional  analysis  of  combined  mode  I - mode  II 
crack.  Special  interpolation  functions  for  element  displace- 
ments are  chosen  so  that  the  displacement  in  planes  normal  to 
the  front  depends  on  the  square  root  of  distance  from  the 
front  and  this  results  in  the  inverse  square  root  stress  mode 
in  these  planes.  Since  the  angular  variations  of  displace- 
ments can  not  follow  those  of  the  corresponding  asymptotic 
terms,  twelve  special  elements  encircle  each  segment  of  the 
crack  front.  For  the  analysis  of  a compact  tension  specimen, 
the  finite  element  mesh  representation  contained  522  elements 
and  1980  degrees  of  freedom.  Stress  intensity  factors  were 
calculated  by  C.O.D.  at  nodes  closest  to  the  crack  front.  He 
gaged  the  candidate  degree  of  error  on  the  stress  intensity 
factor  solution  as  6 per  cent  of  the  exact  solution. 

Blackburn  and  Hellen  [72]  presented  a fifteen-node  penta- 
hedron element  in  which  the  displacement  includes  terms  vary- 
ing linearly  and  as  the  square  root  of  the  distance  from  the 
crack  front.  The  stress  intensity  factors  were  determined  by 
C.O.D.  and/or  by  the  method  of  virtual  crack  extensions  which 
is  similar  to  the  technique  developed  by  Parks.  With  meshes 
containing  between  50  and  100  quadratic  isoparametric  elements 
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and  with  5 minutes  of  CPU  time  on  an  IBM  370/165  computer, 
results  accurate  to  within  a range  of  1 per  cent  to  4 per 
cent  of  known  solutions  were  obtained. 

Raju  and  Newman  [75,76]  used  eight-node  hexahedron  square 
root  elements  which  contain  the  square  root  terms  in  their 
assumed  displacement  field  in  addition  to  the  pentahedron 
singularity  elements  which  are  similar  to  those  of  Tracey's. 

In  their  approach,  8 singularity  elements  encircle  each 
segment  of  the  crack  front,  while  32  square- root  elements 
surround  these  singular  elements.  Their  typical  finite 
element  model  for  the  analysis  of  the  surface  flaw  contains 
4317  degrees  of  freedom.  Mode  I stress  intensity  factor  was 
extrapolated  from  SIF  determined  by  stress  values  at  several 
reference  points  away  from  the  crack  front. 

It  is  well  known  that  an  inverse  square  root  r singular- 
ity occurs  in  isoparametric  finite  elements  if  the  mid-side 
nodes  are  displaced  to  the  quarter  point  from  their  usual 
position  of  any  side.  Barsoum  [67]  extended  the  idea  of  the 
distorted  isoparametric  element  to  three-dimensional  crack 
analysis.  He  constructed  20- node  pentahedron  and  20-node 
hexahedron  special  elements.  However,  later  Hibbit  [71] 
pointed  out  that  the  strain  energy  in  the  distorted  hexahedron 
element  is  unbounded  and  he  recommended  to  use  only  the  penta- 
hedron element. 

Bloom  [69]  studied  the  convergence  of  twenty- node 
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distorted  hexahedron  elements  by  analyzing  a compact  tension 
specimen.  Three  mesh  representations  of  1200,  3528  and  4965 
degrees  of  freedom  were  used  to  model  a quadrant  of  the  speci- 
men. The  resulting  stress  intensity  factors  at  the  midsurface 
were  respectively  -3.5,  6.3  and  7.3  per  cent  higher  than  the 
accepted  two-dimensional  plane  strain  solution.  Both  Barsoum 
and  Bloom  utilized  the  C.O.D.  method  to  obtain  stress  inten- 
sity factor  solutions. 

Recently  Atluri  and  Kathiresan  [66,120,121]  developed  a 
three-dimensional  crack  element  based  on  the  assumed  displace- 
ment hybrid  model.  The  geometry  of  the  crack  element  is  a 
twenty-node  isoparametric  hexahedron  element  with  six  quadra- 
tically  curved  surfaces.  Since  the  stress  intensity  factor 
KI'  KII  an<^  KIII  were  assumed  to  be  quadratic  functions  of 
the  crack  front  coordinates,  the  strain  energy  corresponding 
to  the  singular  strains  had  to  be  evaluated  by  volume  inte- 
grals with  strong  1/r  and  l//r  singularities.  In  their 
initial  attempt  [66]  this  integration  was  done  by  product 
Gaussian  quadrature  using  as  many  as  340  integral  points. 
However,  in  their  later  development  [120]  refined  integration 
schemes  were  employed  to  remove  the  singularities  in  the 
volume  integrals.  Also,  it  appears  that  a curved  crack  front 
element  will  require  the  two  profile  faces  to  be  flat  and 
normal  to  the  curved  crack  front.  In  that  case  there  will, 
in  general,  be  gaps  or  overlaps  in  the  resulting  finite 
element  mesh.  It  is  realized,  however,  that  errors  due  to 


such  finite  element  approximations  will  be  reduced  when  the 
element  size  is  reduced. 

2.5.4  Superposition  of  Analytical  and  Finite  Element 
Solutions 

Yamamoto  and  Tokuda  [20]  proposed  a method  for  the 
stress  analysis  of  two-dimensional  cracks  on  the  basis  of  the 
superposition  of  analytical  and  finite  element  solutions. 
Later,  Yamamoto  and  Sumi  [60]  extended  the  superposition 
method  to  three-dimensional  crack  problems.  Basically,  the 
method  consists  of  arranging  three  solutions,  i.e.  an 
analytical  solution  for  an  infinite  domain,  a finite  element 
solution  which,  when  added  to  the  analytical  solution,  satis- 
fies the  prescribed  external  boundary  condition,  and  a finite 
element  solution  for  a given  boundary  condition  so  as  to  meet 
the  internal  boundary  condition  (stress  free  condition  on 
some  collocation  points  on  the  crack  surface) . The  super- 
position method  is  similar  to  the  Schwartz  alternating  tech- 
nique commonly  used  in  solving  two-dimensional  boundary  value 
problems  for  a doubly  connected  region.  Instead  of  an 
analytical  solution  for  semi- inf inite  domain  required  in  the 
alternating  technique,  a finite  element  solution  is  utilized. 
Therefore,  this  method  has  the  same  drawbacks  as  those  of  the 
alternating  technique.  The  result  is  very  sensitive  to  the 
location  of  the  collocation  points  where  the  stress-free 
conditions  are  enforced  as  well  as  the  analytical  solution 


chosen.  Since  it  requires  an  analytical  solution,  it  can  be 
applied  only  to  cracks  of  simple  geometry. 

A survey  of  existing  finite  element  approaches  for  three- 
dimensional  crack  analyses  is  given  in  Table  1. 


SECTION  3 


HYBRID  STRESS  CRACK  FRONT  ELEMENT  FOR 
THREE-DIMENSIONAL  LINEAR  FRACTURE  PROBLEMS 

3.1  Introduction 

There  exist  two  assumed  stress  hybrid  finite  element 
models  for  the  development  of  special  crack  front  elements 
for  linear  fracture  mechanics.  The  basic  schemes  for  these 
models  are  as  follows: 

(1)  Scheme  I is  based  on  the  complete  expansion  of  the 
exact  near  field  solution  each  term  of  which  satis- 
fies both  equilibrium  and  compatibility  conditions 
inside  the  element  and  the  stress  free  condition 
over  the  c rack  surface.  Independently  assumed 
boundary  displacements  are  assumed  to  satisfy  the 
compatibility  condition  -with  neighboring  elements. 

(2)  Scheme  II  is  based  on  assumed  equilibrating  stress 
field  including  the  dominant  singular  term  as  well 
as  regular  polynomial  basis  term  and  on  independent- 
ly assumed  boundary  displacements.  The  dominant 
singular  term  also  satisfies  the  compatibility 
condition  and  independently  assumed  boundary  dis- 
placements are  assumed  to  satisfy  the  interelement 
compatibility  condition  with  neighboring  elements. 

The  application  of  these  two  schemes  to  two-dimensional 
crack  elements  have  been  reported  by  Tong,  Piar  and  Lasry  [19] , 


by  Pian,  Tong  and  Luk  [12],  respectively.  It  appears  that 
for  two-dimensional  problems  the  most  desirable  element  is 
one  formulated  by  the  first  scheme  for  which  only  one  crack 
element  is  needed  and  an  embedded  crack  is  included  in  the 
element.  In  implementing  such  element,  only  boundary  integra- 
tions are  called  for.  Since  the  element  boundary  does  not 
contain  the  crack  tip,  the  integration  does  not  involve  any 
singular  term. 

For  a three-dimensional  crack  front  element,  it  is  not 
possible  to  obtain  the  complete  expansion  of  the  exact  solu- 
tion in  the  vicinity  of  a crack  front.  The  only  available 
singular  solution  is  the  asymptotically  exact  term  which  is 
derived  from  the  well  known  solution  along  the  plane  normal 
to  the  front  of  an  embedded  elliptical  crack.  Hence,  only 
the  second  scheme  described  above  is  applicable  to  develop 
three-dimensional  crack  front  element. 

In  this  section,  three-dimensional  crack  front  elements 
are  developed  based  on  the  second  scheme  of  the  hybrid  stress 
finite  element  method.  The  basic  concern  in  the  present 
development  is  to  construct  a relatively  simple  finite  element 
model  so  that  it  can  be  used  conveniently  to  solve  everyday 
engineering  problems  in  fracture  mechanics. 

3.2  Outline  of  the  Three-Dimensional  Crack  Element 

The  geometrical  shape  of  the  presently  derived  crack 
front  element  is  an  arbitrarily  oriented  hexahedron  with 


straight  edges  shown  in  Figure  3.  In  order  to  conveniently 
identify  element  faces,  we  call  the  face  1485  and  1265  which 
contain  the  crack  front  line  "frontal  face",  the  face  1234 
and  5678  which  intersect  the  front  line  "profile  face"  and 
the  face  4378  and  2376  which  are  remote  from  the  crack  bound- 
ary "remote  face".  Two  remote  faces  can  be  warped  surfaces 
but  other  faces  should  be  flat.  Since  the  edges  are  all  made 
straight  lines,  if  in  a problem  the  crack  front  is  curved  it 
must  be  approximated  by  an  assemblage  of  straight  segments. 
Each  segment  of  the  crack  front,  then  serves  as  the  common 
edge  of  a group  of  crack  front  elements  that  surround  it.  An 
example  arrangement  is  shown  in  Figure  3.  These  elements  are 
classified  as  type  A element  which  contains  the  crack  surface 
and  type  B element  which  does  not  contain  the  crack  surface. 
Elements  labelled  as  A'  and  B'  are  those  with  the  crack  front 
located  at  the  upper  corners  of  the  block  elements. 

Three  types  of  crack  front  elements,  i.e.,  8-node,  12- 
node  and  16-node  elements  are  constructed  (Figure  6) . The 
16-node  element  has  mid-side  nodes  on  all  edges  on  the  pro- 
file face  that  intersects  the  crack  front,  whereas  the  12- 
node  element  has  mid-side  nodes  only  on  edges  which  are  not 
connected  with  the  crack  front.  Figure  6 indicates  typical 
assemblages  of  crack  front  elements.  The  12- node,  20-node 
and  26-node  half  elements  are  for  problems  which  are  symmetric 
about  the  plane  of  the  crack,  while  the  20-node,  36-node  and 


46-node  superelements  are  for  general  asymmetric  problems. 


It  is  seen  that  since  only  the  two-dimensional  behavior 
of  the  singular  term  is  introduced,  the  elements  are  neces- 
sarily narrow  in  the  direction  along  the  crack  front.  The 
stress  intensity  factors  K^,  and  associated  with 

each  crack  front  element  are  assumed  to  be  constant  and  hence 
to  determine  the  distribution  of  K^.,  and  Kjjj  a number  of 

special  elements  should  be  used  along  the  crack  front. 

3.3  Variational  Formulation 

The  modified  complementary  energy  functional  7rmc  can  be 
written  as  follows  [87], 


fmc  (aij ' ui}  Z [/V  2 Cijk£  aij  °k£  dV 

n=l  n 


+/S  Ti  °i  ds 
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In  the  above  expression, 


- Stress  component 


C. - Elastic  compliance  tensor 

X ^ JC  36 


£L  - Boundary  displacement  component 


T^  - Prescribed  surface  traction  on 
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(3.1) 


V - Volume  of  the  nth  element  (n  = 1,2,3,  . . . , N) 
n 

9Vn  - Entire  boundary  of  the  nth  element 
S - Portion  of  the  element  boundary  where  surface 

O 

n traction  T^  is  prescribed 

N - Total  number  of  finite  elements 

dV  - Elemental  volume 

dS  - Elemental  surface 

T.  - Surface  traction  component 


The  independent  field  variables  are  the  equilibrating 
stresses  and  the  boundary  displacements  u^.  In  the  crack 

front  elements  which  are  directly  adjacent  to  the  crack  front, 
the  assumed  stresses  include  not  only  the  simple  polynomial 
expansion  along  the  three  coordinate  axes,  but  also  the 
asymptotically  exact  terms  which  are  derived  from  the  well 
known  classical  near-field  solutions  of  an  embedded  ellipti- 
cal crack  in  an  infinite  domain  and  are  given  in  Equation 
(2.1).  In  the  far-field  elements  where  the  effect  of  stress 
singularity  is  of  little  importance,  only  regular  polynomial 
basis  terms  are  incorporated  in  the  stress  assumption. 

In  a crack  front  element,  the  assumed  stresses  are 
expressed  as  follows: 


ai-  = a (l//r,  9)  + a (x,y,z)  (3.2) 

ij  ij 

where  ag  corresponds  to  the  singular  stress  solution, 
ij 

Equation  (2.1),  and  a contains  only  regular  polynomials. 

rij 

The  subscript  s and  r refer  to  "singular"  and  "regular" 
respectively.  The  singular  stress  term  satisfies  both  the 
homogeneous  equations  of  equilibrium  and  the  stress  free 
boundary  conditions  over  the  crack  surface,  i.e.. 


0 

0 


in  V 

n 

(3.3) 

on  the  crack  surface 


where  vVs  are  the  components  of  the  unit  normal  vector  to 
the  element  surface.  The  regular  stress  terms  are  selected 
so  as  to  satisfy  the  complete  equations  of  equilibrium  with 
prescribed  body  forces,  i.e. 

a + F.  = 0 (3.4 

rij  > 3 1 

where  iv  is  a component  of  the  prescribed  body  forces  and  a 
known  function  of  coordinates.  The  surface  tractions  for 
crack  front  elements  are  also  divided  into  two  parts  and  are 
related  to  the  interior  stresses  by 


(3.5 


The  regular  stress  term  ar 
parts  as 


ij 


is  further  divided  into  two 


+ o. 


Pij 


(3.6 


where  satisfy  the  homogeneous  equilibrium  equations 

ij 

while  a is  a particular  solution  of  the  equations  of 
pij 

equilibrium  with  the  prescribed  body  forces.  The  subscript 
h denotes  "homogeneous"  while  p denotes  "particular". 


Substituting  Equation  (3.2)  into  the  first  term  of  the 
right  hand  side  of  Equation  (3.1),  the  complementary  strain 


energy  in  the  crack  front  element  is  written  as 
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Since  the  singular  stress  field  og  is  the  exact  asymptotic 
solution  of  an  embedded  elliptical  crack,  it  satisfies  compat- 
ibility conditions  as  well  as  homogeneous  equations  of  equi- 
librium. Hence,  the  volume  integral  of  the  first  term  of 
Equation  (3.7)  can  be  replaced  by  the  surface  integral  through 
the  application  of  divergence  theorem. 
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where  u is  an  asymptotically  exact  displacement  field  which 
si 

corresponds  to  the  singular  stress  field  a . The  explicit 

Sij 

form  of  u is  given  in  Equation  (2.2).  The  regular  stress 
i I 

term  n satisfies  the  equations  of  equilibrium  with  pre- 

ij 

scribed  body  forces  and  hence  the  second  term  in  Equation 
(3.7)  can  be  converted  to 
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Upon  substitution  of  Equations  (3.7),  (3.8)  and  (3.9)  into 

Equation  (3.1),  the  hybrid  variational  functional  it  can  be 


written  in  the  form 
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where  N-  is  the  total  number  of  far  field  elements  and  n is 
f c 

the  number  of  crack  front  elements. 


3.4  Matrix  Formulation 

For  each  crack  front  element  V_  ( n = 1,2,  . ..  N ),  the 

n c 

assumed  stresses  are  written  in  the  matrix  form  as: 
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a - P 8_  + P 6 + PS 

-»  ~s~s  ~ ~ ~p~p 


(3.11) 


where  P 8_  is  the  asymptotically  exact  singular  solution  func- 
~ s ~ s 

tion  given  in  Equation  (2.1)  in  terms  of  r and  0,  8 repre- 

5 

sents  the  vector  of  undetermined  stress  intensity  factors  as 


MrZ- 


6g  - (Kj,  Kjj,  KIIX} 


(3.12) 


In  the  present  formulation,  8 is  assumed  to  be  constant 

within  each  crack  front  element.  P is  the  regular  polynomial 

basis  function  of  spatial  coordinates,  B is  unknown  stress 

parameters  and  PB  satisfies  the  homogeneous  equations  of 

equilibrium.  P 8 is  a known  particular  solution  of  the  non- 

homogeneous  equations  of  equilibrium  and  in  case  any  non- zero 

tractions  are  prescribed  on  the  boundary  of  the  element, 

appropriate  terms  should  be  added  to  PpBp. 

The  surface  tractions  on  the  element  boundary  3Vn 

(n  = 1,2,  ...  N ),  which  relates  to  the  assumed  stress  dis- 
c 

tribution,  can  be  expressed  in  terms  of  the  generalized 
stresses  by  equation  of  the  form: 

T * ?s5s  + ? ! + ?P5p  (3-13) 

in  which  R,  R and  R are  obtained  from  P , P and  P , respec- 

~ S -V  p ~ 5 ~ ~ p 

tively.  On  the  other  hand,  in  the  far  field  elements 
Vm(m  =1,2,  ...  Nf)  only  the  polynomial  basis  terms  are 
included  in  the  stress  and  traction  assumption,  i.e., 


? - ? 5 + ?p!p 

t = r e + Rp6p 


(3.14) 


The  asymptotically  exact  displacement  matrix  which  corresponds 

to  the  singular  stress  matrix  P B_  can  be  written  as: 

~ s ~ s 
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where  u is  the  asymptotic  solution  matrix  of  displacements 
~ s 

defined  by 


{u_  , u„  , u_  } 


(3.16) 


Lg  is  the  exact  displacement  solution  function  given  in 
Equation  (2.2). 

The  assumed  displacement  on  the  element  boundary  3Vn 
can  be  expressed  as: 

u * L q (3.17) 

where 

u = {u^  u2,  u3>  (3.18) 

q is  the  generalized  nodal  displacement  vector  and  L is  the 
interpolation  function  matrix.  In  crack  front  elements,  L 
contains  asymptotically  correct  variation  of  displacements  on 
the  faces  that  meet  or  intersect  the  crack  front  line.  For 

the  other  part  of  the  boundary  of  singular  elements  and  for 

the  entire  boundary  of  far-field  elements,  L consists  only  of 
conventional  polynomial  basis  interpolation  functions. 
Inserting  Equations  (3.11),  (3.13),  (3.14),  (3.15)  and  (3.17) 

into  Equation  (3.10)  gives 
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The  stationary  condition  of  the  above  functional  tt  with 

me 

respect  to  variations  of  stress  parameters  8 which  are 
independent  from  one  element  to  the  other  gives 


H.  . 8 + H, 8^  + H,  8 - G q = 0 for  crack  front  elements 

~hh  - ~hp-p  ~hs~s  ~ ? 

(3.21) 
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for  far  field  elements 


(3.22) 


Equations  (3.21)  and  (3.22)  can  be  solved  for  8 as  follows. 

8 = H.  . ^ (G  q - H,  B - H, 8 ) for  crack  front  elements 

- hh  ~ ~ ~hs~s  ~hp~p 

(3.23) 


8 = H.  . (G  q - H, B„)  for  far  field  elements 

~ nn  - ~ ~hp~p 

(3.24) 

Substituting  Equations  (3.23)  and  (3.24)  back  into  the  func- 


tional Equation  (3.19),  we  obtain  the  modified  complementary 
energy  functional  it  expressed  in  terms  of  generalized  nodal 

ITl  L. 

displacements  q and  the  stress  intensity  factor  vector  B0as: 
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In  the  above,  the  matrix  Hhh  is  a symmetric  matrix.  The  stress 


intensity  factor  8 is  not  completely  independent,  but  is 

S 


common  to  a group  of  crack  front  elements  which  have  one  crack 


front  segment  as  a common  edge.  This  group  of  crack  front 


elements  can  be  conveniently  treated  as  a "superelement". 


Element  matrices  K,  K , K , B'  and  Q for  a superelement  are 

~ ~ SS  ~ i S <v  -v 


constructed  by  summing  the  overlapping  terms  of  the  correspond- 


ing element  matrices  for  the  group  of  crack  front  elements. 


Since  the  constant  term  A'  has  no  effect  on  the  variation,  it 


may  be  dropped  from  the  functional  representation.  Hence,  the 


variational  function  in  Equation  (3.25)  is  rewritten  as: 
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where  N is  the  total  number  of  the  superelements  and  the 


matrices  K,  K , K_,  B'  and  Q denotes  their  assembled  element 
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Since  the  stress  intensity  factor  vector  6 is  independent 

S 

for  different  superelements,  the  vanish  of  the  variation  of 
the  functional  in  Equation  (3.27)  with  respect  to  6e  gives 


«ss  *s  + *rs  ? + !'  - 0 


(3.28) 


The  stress  intensity  factor  vector  8 can  now  be  expressed  in 

~ 2 

terms  of  the  generalized  nodal  displacement  vector  q as 


follows. 
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Substituting  Equation  (3.29)  back  into  Equation  (3.27)  and 

neglecting  constant  terms,  the  modified  complementary  energy 

functional  it  can  thus  be  written  as: 
me 
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where 


T -1 

K = K + K K K 
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In  the  above  expression,  all  the  stress  parameters  are  formally 


eliminated  from  the  functional  and  it  contains  only  generalized 
nodal  displacements  q as  unknown  parameters.  Kg  gives  the 


element  stiffness  matrix  of  a crack  front  superelement  which 
is  compatible  with  those  of  both  assumed  displacement  and 
hybrid  stress  elements. 

The  assemblage  process  of  the  element  matrices  to  the 
global  matrices  also  follows  completely  that  of  conventional 
method.  For  the  purpose  of  completeness,  the  assemblage  and 
solution  procedure  is  briefly  sketched. 

The  element  nodal  displacements  q are  not  completely 
independent  of  those  in  other  elements  because  of  the  require- 
ment of  inter element  compatibility.  They  are  expressed  as  a 
function  of  the  independent  generalized  global  displacements. 
Expressing  all  element  nodal  displacements  q in  terms  of  a 
global  set  q*,  the  functional  irmc  yields 

'me  = I ?*T  5*  9*  + 9*V  (3-32) 

where  K*  is  the  global  stiffness  matrix  obtained  by  assembling 
the  element  stiffness  matrices  associated  with  the  same 
generalized  displacements  for  the  entire  domain.  Q*  is  the 
vector  of  global  generalized  nodal  forces. 

The  stationary  condition  6 TTmc  * 0 with  respect  to  the 
variation  of  q*  yields  the  global  equation 

K*  q*  = Q*  (3.33) 

The  above  equation  is  solved  for  the  global  displacements  q*, 
thus  determining  the  generalized  displacements  q for  each 


49 


element.  Then  as  a final  step  in  the  solution  procedure,  the 
stress  intensity  factor  vector  gg  is  calculated  for  each  crack 
front  superelement  by  substituting  the  corresponding  displace- 
ment solutions  into  Equation  (3.29). 

Instead  of  condensing  the  stress  parameter  (stress  inten- 
sity factor  vector)  8_  a priori,  we  may  carry  it  into  the 

~ 5 

global  system  of  equations.  This  approach  leads  to  a dual 

hybrid  method  with  generalized  displacements  q and  stress 

parameters  8_  as  unknowns.  After  assembling  the  element 

matrices  K,  K , K „ and  Q,  the  variation  of  it  (q*,  B*) 

~ ~rs  ~ s s ~ me  ~ ~s 

leads  to  the  final  algebraic  equation  of  the  form, 


‘k* 

K* 

-rs 

■ 

q*  ■ 

Q* 

k*t 

rs 

(0  * 
(0 

. si. 

0 

(3.34) 


The  above  matrix  equation  can  be  solved  directly  for  the 
generalized  nodal  displacements  q*  and  the  stress  intensity 
factors  6|.  Since  the  global  matrix 


K* 


K*T 

~rs 


K* 
~r  s 


K* 

~ss 


is  not  always  positive  definite,  but  positive  semi-definite, 
this  formulation  may  not  be  compatible  with  the  conventional 
general  purpose  finite  element  program  for  the  assumed  dis- 
placement method. 
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In  the  derivation  of  element  matrices,  the  local  Cartesian 
coordinate  system  (x,y, z)  with  its  origin  being  at  one  end  of 
the  crack  front  edge  is  taken  as  shown  in  Figure  7.  The 


y-axis  is  taken  to  coincide  with  the  crack  front  edge,  while 
the  x-axis  is  taken  to  be  normal  to  the  crack  front  edge  so 
that  the  crack  surface  lies  on  the  x-y  plane.  The  z-axis  is 
normal  to  both  the  x-  and  y-axes.  There  is  no  loss  of 
generality  in  defining  the  local  coordinate  system  in  the 
above  manner.  The  coordinates  (r,  9)  are  defined  on  the  plane 
which  is  perpendicular  to  the  crack  front  and  are  related  to 
the  local  Cartesian  coordinates  as: 

x = r cos9 
z = r SinQ 

The  isoparametric  non-dimensional  coordinate  system  (5,n»0  is 
also  utilized  to  specify  a point  within  the  element.  This 
coordinate  system  also  facilitates  the  numerical  integration 
of  element  matrices.  The  relationship  between  the  non-dimen- 
sional coordinate  (£,n,0  at  any  point  P and  the  local 
Cartesian  coordinate  (x,y,z)  can  be  written  as 

8 

x = z <f> . (C,n#  Ox. 

i=l  1 1 

8 

y = E <j> . (5, n, Oy, 

i=l  1 1 


jtf  . 1 ZT  ZZ 


where 

4>i  = | (l+qO  (l+Hj^n)  (1+^5) 

(x.,  y^,  z.)  and  (£^,  are  the  Cartesian  coordinate  and 

the  non-dimensional  coordinate  of  the  ith  corner  node  of  the 
crack  front  element.  The  non-dimensional  coordinate  system 
(£»n»S)  are  scaled  so  that  six  element  faces  are  defined  by 
£ = ±1,  n = ±1,  £ = ±1.  For  example,  for  a type  B element  in 
Figure  3,  the  profile  faces  determine  n = ±1,  the  frontal 
faces  define  £ = -1,  £ = -1  and  the  remote  faces  are  £ = 1, 

C - 1. 

For  the  twelve- node  and  sixteen- node  crack  front  elements, 
whose  geometrical  shape  is  a hexahedron  with  straight  edges, 
the  mid-side  nodes  do  not  affect  the  original  element  shape, 
i.e.,  the  mid-side  nodes  do  not  affect  the  relation  between 
the  non-dimensional  and  the  Cartesian  coordinates.  Therefore, 
for  the  twelve-node  and  sixteen-node  element,  the  same  rela- 
tionship as  that  in  Equation  (3.35)  holds  between  these  two 
coordinate  systems. 

3.6  Assumed  Stress  Field 

The  stresses  to  be  assumed  inside  the  element,  which  are 

a , a , a , t , t and  t , should  satisfy  the  following 
x y z y z zx  xy 

homogeneous  equations  of  equilibrium. 


... 
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Since  the  prescribed  boundary  stresses  do  not  constitute  a 
restrained  boundary  condition  in  the  hybrid  stress  finite 
element  model,  it  is  not  necessary  to  select  the  stress  field 
so  as  to  satisfy  the  natural  boundary  conditions.  However,  in 
the  type  A and  type  A'  crack  front  elements  which  contain  the 
crack  surface,  the  assumed  stresses  are  chosen  to  satisfy  the 
stress  free  condition  over  the  crack  surface  for  obtaining 
more  accurate  stress  intensity  factor  solutions. 

The  magnitude  of  the  stress  intensity  factor  and 

Ki:ri  will  vary  from  point  to  point  along  the  crack  front. 
However,  in  this  case,  it  is  impossible  to  construct  the  self 
equilibrating  singular  stress  field  with  K ^ , K and  KII;I 

being  simple  functions  of  the  crack  front  coordinate  y,  and 
hence,  in  the  present  formulation,  the  stress  intensity 
factors  are  assumed  to  be  constant  within  each  crack  front 
element.  The  singular  part  of  the  assumed  stresses  with 
constant  Ki;[  and  K is  listed  in  Table  10.  Since  it  is 
the  exact  asymptotic  solution  for  a line  crack  embedded  in  an 
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infinite  domain,  it  satisfies  both  the  homogeneous  equations 
of  equilibrium  in  the  element  interior  and  the  stress  free 
condition  over  the  crack  surface. 

Although  no  method  has  been  reported  to  date  to  determine 
the  sufficient  and/or  optimum  number  of  stress  parameters 
a priori,  the  number  of  stress  parameters  cannot  be  chosen 
freely.  The  following  necessary  condition  should  always  be 
met  to  avoid  kinematic  deformation  modes  which  involve  zero 
strain  energy. 

m > n-Jl  (3.37) 


where  m is  the  total  number  of  stress  parameters  6 and  3g, 
n is  the  total  number  of  generalized  nodal  displacements  and 
l is  the  number  of  rigid  body  modes.  For  a general  three- 
dimensional  problem  £ = 6. 

The  regular  polynomial-based  self-equilibrating  stress 
field  is  derived  from  the  three-dimensional  stress  functions. 
First,  Maxwell's  stress  functions  X2  anc^  X3  are 

assumed  in  the  polynomial  form  in  terms  of  Cartesian  coordi- 
nates x,  y and  z.  Then,  all  stress  components  are  computed 
algebraically  based  on  the  following  stress  representations. 


a 


x 


3 


2 

X 


2 


a 


y 
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.2  ~ 2 
3 X2  . 9 xi 

az  " — r*  + ~r 


9x 


3y 


y z 


zx 


xy 


- 2 

3 *1 
9y9  z 

~ 2 

3 x2 

9z  9x 
~ 2 

3 X 3 
9x3y 


(3.38) 


In  type  A and  type  A'  crack  front  elements,  the  stress  field 
derived  in  the  above  mentioned  manner  should  be  modified  so 
as  to  satisfy  the  stress  free  condition  over  the  crack  surface 
(z  = 0)  . 

Within  type  B and  type  B'  8-node  crack  front  element,  the 
stress  field  is  assumed  as 


X 

B1  + 

e7y  + 

BgZ  + 

B9X  + B10 

yz 

y 

62  + 

eux  + 

B12Z 

+ 

Bny  + 

®142x 

z 

63  + 

Bi5x  + 

B16y 

+ 

S17Z  + 

S18xy 

yz 

= B4  - 

(617 

+ B20 

>y 

+ S19z 

+ 822x 

zx 

= B 5 “ 

(Sg  + 

B21> 

z + 820x  + 

®23y 

xy 

= B6  - 

(6U 

+ b19 

) x 

+ S21y 

+ 824x 

(3.39) 
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Within  type  A and  type  A'  8-node  element,  the  stresses 
are  assumed  as 

°x  “ 61  + B4y  + S5Z  + S6X  + S13Xy  + S14yz  + B15ZX 


+ 

( B19 

+ 820)x2  + B22x2y 

+ e24x2z 

II 

>1 

o 

e2  + 

87  Z + BgX  + B9y  + 

B16yz  + B1  i7*  + B18xy 

+ 

(620 

+ S21)y2  + B23xy2 

+ «24y2z 

az  = 

(S19 

+ B21)z2  + B22yz2 

+ 623z  x 

yz 

" (69 

- B10> 

z - 618zx  - 

2B21yz  " 

2B23XyZ 

zx  — 

- (66 

' Bll} 

z - B13yz  - 

2B19zx  - 

2B22xyz 

xy 

B3  “ 

Biox  - 

8lly  + B12Z 

- S15yz  - 

B16ZX  " 2B20Xy 

- 2824xyz  (3.40) 

The  above  stress  field  satisfies  the  stress  free  condition 
over  the  crack  surface  as  well  as  the  homogeneous  equations 
of  equilibrium. 

Regular  polynomial-based  assumed  stress  field  used  in  the 
present  study  are  listed  in  matrix  form  in  Tables  2,  4,  6 and 
8.  Table  14  also  lists  the  number  of  B's  used  in  the  various 
crack  front  elements. 


n 

I 

1 4, 
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3.7  Body  Forces 

Body  forces  may  be  caused  by  gravity,  acceleration, 
magnetic  field  and  so  on.  The  equivalent  nodal  force  vector 
associated  with  the  prescribed  body  forces  is  obtained  as 


follows , 


For  a far  field  element 


9b  - * Sp§p  + ST  Shh'1  5hP!p 


(3.41) 


For  a crack  front  element 


rp  iT»  — 1 rn  _ 1 rp  —1 

Q = - G B + G H,  , H.  3_  “ K"1  K (B  - H7  H.  , H,  6 ) 

~B  ~p~p  - ~hh  ~hp~p  ~rs  ~ss  ~ ~hs  ~hh  ~hp~p 

(3.42) 


To  insure  a unique  finite  element  solution,  there  is  a require- 
ment for  the  given  particular  solution  that  the  order  of  the 
polynomials  in  the  particular  solution  should  be  equal  to  or 
lower  than  that  of  the  highest  complete  polynomials  in  the 
homogeneous  stress  field  [87].  For  type  B and  B'  crack 
elements  and  far-field  elements,  the  following  particular 
solution  can  be  used  for  constant  body  forces 


a = 
x 


F x 
x 


ay  = 


Fy  y 


a2  = 


F z 
z 


T = T 

zx  xy 


(3.43) 


However,  for  type  A and  A'  crack  front  elements,  if  the  stress 

assumption  given  in  Equation  (3.40)  is  used,  the  assumed  stress 

a , a and  a lacks  the  constant  term  and  hence  a unique 
z yz  zx 

solution  is  not  guaranteed  for  the  particular  solution  of  any 
order.  In  this  case,  the  equivalent  nodal  forces  may  be  cal- 
culated in  the  inconsistent  manner. 

3.8  Assumed  Boundary  Displacement  Field 

A crack  front  element  is  connected  to  conventional  solid 
elements  on  its  remote  faces,  while  it  joins  to  other  singular 
elements  through  the  frontal  and  profile  faces.  For  a success- 
ful implementation  of  this  crack  front  element,  displacements 
must  be  compatible  along  those  interelement  boundaries.  This 
condition  is  normally  satisfied  in  the  hybrid  model  by  assum- 
ing boundary  displacements  in  terms  of  interpolation  functions 
and  generalized  nodal  displacements  on  each  element  interface. 
For  far-field  (regular)  elements,  a single  interpolation  func- 
tion can  be  used  for  the  assumptions  of  all  displacement  com- 
ponents over  all  element  faces  if  each  face  consists  of  the 
same  number  of  nodes  and  same  geometrical  shape.  However,  in 
a crack  front  element,  one  has  to  work  with  a different  dis- 
placement interpolation  for  each  face  to  simulate  the  correct 
asymptotic  behavior  of  displacements  in  the  neighborhood  of 
the  crack  front.  In  the  present  study  the  boundary  displace- 
ments adopted  follow  the  approach  used  by  Atluri  and 
Kathiresan  [66,120],  in  connection  with  their  crack  element 
by  the  assumed  displacement  hybrid  model. 


Assumed  Boundary  Displacements  for  8-Node  Element 


Figure  7 shows  the  typical  8-node  crack  front  element 
and  its  relative  node  numbering  scheme. 

Frontal  Face  (1,  4,  8,  5,  - 1,  2,  6,  5) 

On  the  frontal  face,  which  contains  the  crack  front  edge, 
the  asymptotic  behavior  of  displacements  in  Equation  (2.2) 
indicates  square  root  r variation.  Displacement  component 
u,  v,  w can  use  the  same  interpolation  polynomials.  For 
example,  on  face  1,  4,  8,  5 (C  = -1),  the  x-direction  dis- 
placement u is  assumed  as 

5 = 3^  + a2t/r  + a3n  + a4/r  n 

where  all  a's  are  constants. 

Prescribing  the  nodal  displacement  in  terms  of  the  nodal 
values  of  coordinates  r and  n,  one  obtains  a set  of  simulta- 
neous equations  as 

+ a2/rT  + a3ni  + a4/r7  (i  = 1,4, 8, 5) 

where  r^,  are  the  coordinates  r and  n of  the  ith  node. 
Solving  these  equations,  the  parameters  a^(i  = 1,2, 3, 4)  are 
uniquely  expressed  in  terms  of  the  corresponding  generalized 
nodal  displacements,  viz,  q^,  q1Q,  <322'  q13‘  Ttie  ^ and  z~ 
direction  displacement  v and  w are  assumed  in  a similar 


fashion. 


Profile  Face  (1,  2,  3,  4,  - 5,  6,  7,  8) 

On  the  profile  face  which  intersects  the  crack  front 
edge,  the  asymptotic  behavior  of  displacements  takes  the  form 
of  a trigonometric  function  of  9 bounded  by  square  root  r as 
shown  in  Equation  (2.2).  To  simulate  this  correct  angular 
and  radial  variation  of  displacement,  the  boundary  displace- 
ments over  this  face  are  assumed  as 

0 30 

u = a^  + a2»/r  [ ( 2k— 1)  cos  2 - cos  — ] 

0 30 

+ a 3 /r  [(2k+3)  sin  2 + sin  y-]  + a4r 

0 0 

v = b^  + b2/r  sin  + b3/r  cos  + b4r 

0 30 

w ■ c^  + c2/F  [(2k+1)  sin  - sin  ] 

0 30 

+ c^/r  [(2<-3)  cos  2 + cos  j-]  + c4r 

where  the  coordinates  r and  0 are  measured  on  the  plane 
perpendicular  to  the  crack  front  line,  k is  related  to  the 
Poisson's  ratio  by  < = 3-4v. 

Inserting  nodal  values  of  r and  9 into  the  above  equa- 
tions and  equating  them  to  appropriate  generalized  nodal 
displacements,  we  obtain  three  sets  of  simultaneous  equations. 
Upon  solving  these  equations,  all  coefficients  (a^,a2,a3,a4) , 
(b1,b2,b3,b4)  and  (c1,c2,c3,c4)  are  evaluated  uniquely  in 
terms  of  their  respective  generalized  nodal  displacements. 
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For  example,  on  face  1,  2,  3,  4 (a^a^a^a^  are  determined 
in  terms  of  (q^^, q4, q?/ q1Q) , (b,  ,b,,b7,b4)  are  determined  by 


(q2,q5,q8,q11)  , and  (c^c^CyC^)  are  determined  by  (q3,qgfq 


q12}  ‘ 

When  the  crack  front  is  straight,  boundary  displacement 
is  compatible  with  that  of  the  adjoining  crack  front  element 
on  this  face.  However,  when  the  crack  front  is  curved,  the 
hypothetically  rectilinear  crack  front  line  should  be  symmet 
rical  across  the  profile  face  to  maintain  the  interelement 
displacement  compatibility  on  this  face. 

Remote  Face  (2,  3,  7,  6,  - 4,  3,  7,  8) 


On  the  remote  face  which  is  distant  from  the  crack  front 
the  well  known  interpolation  function  of  polynomial  form  can 
be  used  for  boundary  displacements  assumption. 

For  example  on  face  2,  3,  7,  6,  boundary  displacements 
are  interpolated  as 


u * <(>1  q4  + <t>2  <3 7 + $3  319  + $4  ^16 


v = <^1  q5  + <D2  q8  + 3 q20  + *4  q17 


W - 4>x  q6  + <f>2  q9  + <J>3  q21  + <f>4  q18 


where 


<j>l  = + d-5)  (l-n) 


<t,  = t (1+5)  d-n) 


r 


<f3  = -J  (1+C)  (l+n) 


4>4  = I (1-;)  (l+n) 


It  is  obvious  that  when  a 8-node  far-field  element  based  on 
the  assumed  displacement  model  is  connected  at  this  face, 
interelement  compatibility  of  displacements  is  preserved. 
However,  strictly  speaking,  the  assumed  boundary  displacement 
is  discontinuous  across  the  element  edge  2-3,  3-4,  6-7  and 
7-8.  A number  of  performance  tests  which  will  be  described 
later  shows  that  this  slight  inconsistency  has  no  appreciable 
effect  on  the  finite  element  solution. 

Assumed  Boundary  Displacements  for  12-Node  Element 
A typical  12-node  crack  front  element  and  its  relative 
node  numbering  scheme  are  shown  in  Figure  8. 

Frontal  Face  (1,  4,  8,  5,  - 1,  2,  6,  5) 

Since  there  is  no  intermediate  side  node  on  this  face, 
the  same  interpolates  are  used  as  those  on  the  frontal  face 
of  8-node  element. 

Profile  Face  (1,  2,  9,  3,  10,  4,  - 5,  6,  11,  7,  12,  8) 
Since  there  are  two  additional  mid-side  nodes  on  this 
face,  displacements  are  assumed  as 

Q O 0 

u = a1  + a2>/r  [(2ic-l)  cos  y - cos  y—  ] 


0 30 

+ a3/r  [(2<+l)  sin  y - sin  y-H 


i 
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0 30 

+ ag/r  [ (2k-3)  cos  + cos  j-]  + agr 


v and  w are  given  by  identical  expression  with  (a^,a2,  ...  ag) 
above,  being  replaced  by  (b^,b2,  ...  bg)  and  [c^,c2,  •••  cg)  , 


respectively. 

Remote  Face  (4,  10,  3,  7,  12,  8,  - 2,  9,  3,  7,  11,  6) 
For  example  on  face  4,  10,  3,  7,  12,  8,  displacements 
are  interpolated  as 


Q " 4*1  q10  + ^2  q7  + ^3  q19  + ^ 4 q22  + ^5  q28  + ^6  q34 

V = ^i  qn  M2  q8  + <t> 3 q20  + ^4  q23  + ^ 5 q29  + *6  q35 

W = 4>1  q12  + <P2  q9  +.4>3  q2i  + ^4  q24  + S q30  + *6  q36 

where 

<P1  = - J d-5)  d-n)C 

<p2  = \ (1+0  (l-nK 

<P3  = \ (l+O  (1+nK 

i 

4>4  = ~ \ (1-5)  (1+nK 

<t>5  = \ (l-C2)  (l-n) 


In  the  present  formulation,  the  evaluation  of  element 
matrices  entails  the  integration  of  both  surface  integrals 
and  volume  integrals.  Some  of  these  surface  integrals  contain 
singularities  at  the  crack  front.  This  section  gives  a brief 
comment  on  the  integration  techniques  utilized  in  the  present 
study • 

Any  face  of  the  crack  front  element  is  conveniently 
represented  by  a set  of  three  parametric  equations.  For 
example,  the  element  face  n = 1 in  Figure  7 is  expressed  in 
the  form  of 

x = x (£, ;) 

y = y (£»£) 

z = z (5,0 

The  differential  surface  dS  on  this  face  is  expressed  in  terms 
of  d£  and  d?  as 


where 


r = xi  + yj  + zk 


Expanding  the  vector  cross  product  leads  to 

✓a(C,0  = (a  l + a2  + a*) 1/2 

where 
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laU  5)J 


n=i 
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3 (z,x)  . 
3(C,5)  J 


n=i 


(x,y) 
, C) 


n=l 


3y  3z  3y  3z 
l3£  35  " 35  35*, 


3z  3x  3z  3x 
3£  3C  " 3C  9C  n=1 


3x  3y  3x  3 y 
l3C  3 5 ” 35  3?  n=1 


Thus,  the  surface  integral  of  an  arbitrary  function  f(£,5) 
over  the  face  n = 1 can  be  written  as 

//  x f(5,C)dS  = f(5»C)  /aTcTTT  dCd5 

This  idea  of  evaluating  surface  integral  is  similar  to  the 

idea  used  by  Atluri  et  al.  [66,120]. 

In  the  evaluation  of  the  strain  energy  due  to  the  exact 

asymptotic  stress  and  displacement  field,  the  matrix  Rs  in 

Equation  (3.13)  contains  a 1/ /r  variation  over  the  surface 

while  L contains  a /r  variation.  Therefore,  the  matrix  H 

~ s ~ ss 

in  Equation  (3.20)  is  completely  free  from  singularity  and 
the  aforementioned  surface  integration  technique  can  immedi- 
ately be  applied  to  its  numerical  integration.  On  the  other 
hand,  the  integrand  of  Gs  in  Equation  (3.20)  still  involves 
a l//r  singularity  at  the  crack  front.  When  Gs  is  evaluated 
on  both  the  frontal  and  profile  faces  which  include  the  crack 


front  and  the  singularities,  proper  changes  of  the  integration 
variables  are  necessary  in  order  to  improve  the  numerical 
accuracy. 

On  the  frontal  face,  for  example,  on  the  face  5 = -1, 
only  the  variable  £ is  transformed  into  a new  variable  t by 

C = \ (t+1) 2 - 1 

The  surface  integral  can  be  rewritten  as 

f(n,c)  /a(n,?)  dnd? 

= f(n,s(t))  /a  (ri,  5 (t) ) (t+1)  dndt 

In  the  numerical  integration  of  Gg  over  this  frontal  face, 
3-point  Gaussian  quadrature  was  used  for  n while  5-point 
quadrature  was  used  for  the  new  variable  t. 

On  the  profile  face  n = +1,  both  5 and  £ are  transformed 
into  new  variables  s and  t by 

C = § (s+1) 2 - 1 

C = \ (t+1)2  - 1 

The  surface  integral  over  this  face  is  written  as 

l 

/-1/-1  f(S'^  /a(s,o  dCd?  ' 

= f (?(s),5(t))  /a(5(s)  (tyy  (s+1)  (t+1)  dsdt 


On  this  face,  numerical  integration  was  performed  in  terms  of 
s and  t by  using  5x5  point  Gaussian  quadrature. 

On  the  remote  face,  any  integrand  does  not  contain 
singularity  and  hence  3x3  point  Gaussian  quadrature  was  used 
for  the  evaluation  of  surface  integrals  over  this  face. 

The  volume  integrals  in  the  present  model  are  completely 
free  from  singularities.  The  differential  volume  dv  is 
expressed  in  terms  of  d£,  dn  and  d£  as 

dV  = |jlabs  d?dnd? 

where  |j|  is  the  Jacobian  defined  by 
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The  volume  integral  is  written  as 

fff  f (x, y, z)  dV 
n 

= f (xU,n,c),y(e,n»0rz(5,Ti,0)  |J|dCdnd^ 

For  the  evaluation  of  the  matrix  in  Equation  (3.20),  one 

can  reduce  the  computation  time  in  the  following  manner. 
Every  element  of  Shh'  which  consists  of  an  integral  of  a 
linear  combination  of  a small  number  of  distinct  terms  such 
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as  x,y,z,x  , y ,z  ,xy  and  so  on,  is  expressed  algebraically 
a priori.  Numerical  integration  needs  be  carried  out  only 
for  these  distinct  terms.  For  example,  in  an  8-node  type  B 
element,  only  159  entires  of  the  upper  or  lower  triangular 
part  of  Hhh  are  not  identically  zero.  In  these  nonzero 

elements,  there  are  only  23  distinct  terms,  i.e.,  l,x,y,z,xy, 

222  22  22  22  222222 
yz,zx,x  ,y  ,z  ,xy  ,x  y,yz  ,y  z,zx  ,z  x,xyz,x  y ,y  z ,z  x , 

2 2 2 

xyz  , xy  s,x  yz.  In  an  8-node  type  A element,  the  upper  or 
lower  triangular  matrix  of  has  244  nonzero  entries,  which 

consists  of  a linear  combination  of  72  distinct  terms.  These 

2 2 2 2 2 2 
terms  are  l,x,y, z,  xy,yz, zx,x  ,y  ,z  ,xy  ,x  y,yz  , 

2 223332  2 233  33  33 

y z,zx  ,z  x,x  ,y  ,z  ,x  yz,xy  z,xyz  ,xy  ,x  y,yz  ,y  z,zx  ,z  x, 
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x y ,y  z ,z  x ,x  ,y  ,z  ,x  y,xy  ,y  z,yz  ,z  x,zx  ,xyz  ,x  yz, 

3 222222  2332322332234 

xy  z,xy  z ,x  yz  ,x  y z,x  y ,x  y ,y  z ,y  z ,z  x ,z  x ,x  yz, 

4 4 23  3223  32  233222242 

xy  z,xyz  ,xy  z ,xy  z ,xy  z,xy  z,xyz  ,xyz  ,xy  z ,xy  , 
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x Y ' y z , y z , z x , z x , x y , y z , z x . 


3.10  On  the  Crack  Front  Element  with  Curved  Boundaries 

One  of  the  most  desirable  geometrical  shape  for  three- 
dimensional  crack  front  elements  will  be  a general  hexhahedron 
with  arbitrarily  curved  boundaries.  Figure  4 shows  such  a 
general  "brick"  element  located  at  the  crack  front.  A certain 
group  of  practical  brick  elements  with  curved  boundaries  can 
easily  be  constructed  by  using  the  well-known  isoparametric 
transformation  technique.  However,  because  of  the  necessity 
of  incorporating  the  two-dimensional  asymptotic  behavior 


in  the  finite  element  formulation,  the  actual  shapes  of  the 
hexagonal  elements  must  be  somewhat  restricted.  The  only 
available  asymptotic  solution  near  a crack  front  is  a two- 
dimensional  behavior  in  terms  of  r and  9 , where  r and  9 are 
the  polar  coordinates  with  the  crack  tip  as  their  origin. 

Let  A and  B be  the  two  arbitrary  points  on  the  profile 
face  which  intersects  the  crack  front  line.  If  the  face  is 
warped,  the  two  points  are,  in  general,  on  two  different 
normal  planes  perpendicular  to  the  crack  front  line.  This 
indicates  that  there  is  no  one  to  one  (unique)  correspondence 
between  the  two-dimensional  coordinates  (r,9)  and  an  arbitrary 
point  on  the  warped  profile  face.  Hence,  on  these  element 
faces,  it  is  impossible  to  construct  a reasonable  boundary 
displacement  interpolation  function  which  contains  asymptot- 
ically correct  variation  of  displacement.  Therefore,  profile 
faces  are  required  to  be  flat.  In  addition  to  that,  if  the 
front  edge  of  the  crack  element  is  curved  and  the  profile 
face  is  flat  but  not  normal  to  the  front  line,  the  normal  to 
the  crack  front  which  passes  through  an  arbitrary  point  A on 
the  profile  face  may  intersect  the  crack  front  line  outside 
the  element.  Asymptotic  behavior  of  displacements  and 
stresses  at  point  A can  be  expressed  only  in  terms  of  the 
coordinates  r and  9 and  their  origin  lies  on  the  curved  front 
outside  the  element.  It  is  not  practical  to  calculate 
element  matrices  by  using  the  structural  geometries  outside 
the  element.  Hence,  if  the  crack  front  edge  is  curved,  the 
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element  profile  faces  should  be  flat  and  normal  to  the  cracjc 
front.  However,  there  arises  another  inconsistency.  The 
higher  order  elements  with  curved  boundaries  are,  in  general, 
derived  by  the  nonlinear  coordinate  transformation.  For 
twenty-node  general  brick  element,  the  edge  line  is  made  of 
quasi-quadratic  (parabolic)  curve.  Thus,  if  we  analyze  a 
curved  crack,  for  example  a penny  shape  crack,  by  using 
twenty- node  and/or  other  higher  order  elements,  the  circular 
crack  front  is  approximated  by  the  assemblage  of  parabolic 
and/or  other  incomplete  polynomial  basis  segments  which  do 
not  provide  a complete  fit  to  the  actual  circular  crack  front. 
However,  the  profile  face  of  such  elements  should  be  flat  and 
normal  to  the  element  front  edge,  so  the  overlapping  or 
gapping  is  unavoidable  between  two  neighboring  crack  front 
elements  (Figure  5) . It  is  realized  that  errors  due  to  such 
approximations  in  finite  element  formulations  will  become 
small  when  the  element  size  becomes  small.  The  main  reason 
for  the  choice  of  element  with  a straight  crack  front  is  the 
ease  of  satisfying  the  equilibrium  conditions  for  the  singular 
stress  terms. 

3.11  Evaluation  of  Crack  Front  Elements 

To  assess  the  accuracy  of  the  crack  front  elements 
developed  in  the  present  study,  a simple  test,  following 
that  given  by  Refs.  66  and  120  has  been  performed.  In 
this  and  all  the  following  numerical  computations,  the 
IBM  370/165  at  the  Information  Processing  Center  of  M.I.T. 
has  been  used.  For  each  crack  front  element,  the  nodal  dis- 
placements are  calculated  from  the  known  analytical  solution 


viz..  Equation  (2.2)  that  corresponds  to  the  mixed  mode 
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condition  of  Kj  = KI;J.  = KII3.  = 1.00.  Using  these  displace- 
ments as  input  data,  the  values  of  Kj , K.^  and  are  then 

computed  from  the  numerically  generated  element  matrices  by 
using  Equation  (3.29).  The  results  obtained  for  different 
types  of  crack  front  elements  and  for  different  half  elements 
and  superelements  are  listed  in  Table  14.  It  is  seen  that 
they  all  compare  favorably  with  the  correct  values.  Obviously 
elements  with  more  nodes  yield  more  accurate  solutions. 

3.12  Numerical  Examples  of  Three-Dimensional  Cracks 

3.12.1  Rectangular  Prism  with  a Single  Edge  Crack 
Figure  10  shows  a rectangular  prism  of  width  b thickness 
2t  and  length  2L  containing  a single  edge  crack  of  depth  a. 

A uniformly  distributed  normal  traction  is  applied  over  the 
upper  and  lower  surfaces  which  are  parallel  to  the  crack 
plane.  The  orientation  of  the  global  coordinate  axes  are 
drawn  in  Figure  10  with  the  origin  at  the  center  of  the  crack 
front;  the  Y-axis  is  taken  along  the  crack  front,  the  X-Y 
plane  lies  on  the  crack  surface,  and  the  Z-axis  is  normal  to 
the  crack  surface.  Theoretical  value  of  the  stress  intensity 
factor  Kj  for  a single  edge  crack  in  a rectangular  plate  has 
been  obtained  by  Bowie  and  Gross  based  on  two-dimensional 
plane  strain  assumption  and  is  found  in  Reference  33. 

Since  the  geometry  and  the  loading  is  symmetric  across 


the  plane  Y=0  and  Z=0,  only  a quadrant  of  the  prism  is 
required  for  finite  element  analysis.  Figure  11  shows  two 


different  finite  element  mesh  representations.  In  both  mesh 
representations,  a quadrant  of  the  prism  is  divided  into  five 
equal  layers  by  planes  perpendicular  to  the  crack  front.  Each 
layer  contains  one  upper  half  superelement . The  mesh  in 
Figure  11(1)  is  constructed  from  30  regular  elements  and  5 
super elements  or  a total  of  35  elements.  The  nodes  are 
arranged  at  the  following  locations. 

X/b  = -0.5,  -0.0625,  0,  0.0625,  0.5 

Z/L  = 0,  0.0125,  1 

The  mesh  in  Figure  11(2)  is  constructed  from  80  regular 
elements  and  5 super elements  or  a total  of  85  elements.  The 
nodes  are  placed  at  the  following  locations. 

X/b  = -0.5,  -0.28125,  -0.0625,  0,  0.0625,  0.28125,  0.5 

Z/L  = 0,  0.125,  0.5625,  1 

12-node  upper  half  superelement  is  used  for  two  different 
finite  element  meshes,  while  20- node  superelement  is  used  only 
for  the  coarse  mesh  model.  For  the  remaining  elements,  8-node 
hexahedron  elements  based  on  the  hybrid  stress  model  are  used. 
However,  when  20-node  superelements  are  used,  special  10-node 
hybrid  stress  hexahedron  elements  encircle  each  superelement 
to  connect  8-node  far-field  elements  to  20-node  superelements. 
When  12-node  superelements  are  used,  the  finite  element  model 
of  coarse  mesh  contains  90  nodes,  270  degrees  of  freedom, 
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while  it  contains  114  nodes,  342  degrees  of  freedom  when  20- 
node  superelements  are  used.  In  the  fine  mesh  subdivision 
there  are  168  nodes,  504  degrees  of  freedom.  To  simulate  the 
effects  of  the  adhacent  prism  quadrants,  symmetric  boundary 
conditions  are  imposed  on  the  plane  Y=0  and  Z=0.  The 
Poisson's  ratio  v = 0.3  is  used  in  this  case. 

The  numerical  values  of  the  stress  intensity  factor  Kj 
thus  obtained  are  listed  in  Table  15.  The  variation  of  KI 
solution  along  the  crack  front  and  its  comparison  with  the 
two-dimensional  plane  strain  solution  are  shown  in  Figure  12. 

The  values  of  K^.  in  Table  15  and  in  Figure  12  are  normalized 

. 1/2 
by  the  exact  two-dimensional  plane  strain  solution  K^.  = aa 

for  a crack  of  length  2a  in  an  infinite  plate.  The  KT  solu- 

tion  takes  the  highest  value  in  the  superelement  closest  to 

the  center  of  the  prism  and  gradually  decreases  as  the  surface 

is  approached.  Two-dimensional  plane  strain  solutions  obtained 

by  Bowie  and  Gross  for  crack  length  a,  total  width  b and  total 

1/2 

length  2L  are  expressed  in  the  form  of  Kj  = a (air)  ' F(a/b), 

where  a is  the  applied  stress.  For  a/b  = 0.5,  L/6  = 1,  the 
numerical  values  of  F(a/b)  obtained  by  Bowie  and  Gross  are 
2.86  and  2.82,  respectively.  When  the  fine  (85  elements) 
esh  subdivision  is  used,  in  the  superelement  adjacent  to  the 
center  plane,  the  12-node  superelement  gives  a value  of 
F(a/b)  of  2.89  which  is  1%  higher  than  Bowie's  plane  strain 
solution  and  is  2.5%  higher  than  Gross's  solution.  When  the 


20-node  superelement  is  used,  the  coarse  (35  elements)  mesh 
representation  yields  a value  of  F(a/b)  of  2.76  which  is  3.5% 
lower  than  that  of  Bowie  and  2.1%  lower  than  that  of  Gross. 

It  is  seen  that  the  value  of  F(a/b)  obtained  by  the  combina- 
tion of  12- node  superelement  and  the  coarse  mesh  subdivision 
is  still  within  the  needed  engineering  accuracy. 

3.12.2  Rectangular  Prism  with  a Center  Crack 

Figure  10  shows  a rectangular  prism  of  width  2b,  thick- 
ness 2t  and  length  2L  containing  a center  crack  of  length  2a. 
The  origin  of  the  global  Cartesian  coordinates  (X,Y,Z)  is 
taken  at  the  center  of  the  prism.  A uniformly  distributed 
tension  load  is  applied  on  both  upper  and  lower  surfaces  of 
the  prism.  A two-dimensional  plane  strain  solution  of  the 
stress  intensity  factor  K^.  for  a center  crack  in  a rectangular 
plate  has  been  obtained  by  Ishida,  by  extensively  developing 
the  mapping  functions  and  his  results  are  found  in  Reference 
33. 

Since  the  geometry  and  the  loading  are  symmetric  across 
the  plane  X=0,  Y=0  and  Z=0,  only  one  eighth  of  the  prism  needs 
be  modelled  for  finite  element  analysis.  The  same  mesh  repre- 
sentations as  those  in  the  preceding  example  are  again  used 
in  this  problem.  Symmetric  boundary  conditions  are  imposed 
on  the  plane  X=0,  Y=0  and  Z=0. 

The  stress  intensity  factor  K ^ solution  thus  obtained 

1/2 

are  shown  in  the  normalized  form  Y = KI/ca  ' 


in  Table  16  and 
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in  Figure  13.  Ishida's  two-dimensional  solutions  which  are 

considered  to  be  accurate  within  1 or  2 per  cent  are  expressed 

1/2 

in  the  form  of  KI  = a (air)  F(a/b),  where  a is  the  applied 
stress,  a is  the  half  crack  length,  b is  the  half  plate  width. 
The  a/b  ratios  used  by  Ishida  do  not  exactly  match  that  of 
the  present  example.  His  F (a/b)  values  are  1.18  for 
a/b  = 0.466,  1.25  for  a/b  = 0.535  and  1.33  for  a/b  = 0.592. 

When  the  fine  (85  elements)  mesh  representation  is  used,  the 
use  of  the  12-node  superelements  gives  a F(a/b)  value  of  1.18 
in  the  element  adjacent  to  the  center  plane  for  a/b  = 0.5. 

In  case  of  the  coarse  (35  elements)  mesh  representation,  the 
F(a/b)  value  obtained  by  using  the  20-node  superelement  is 
1.21.  These  F(a/b)  values,  correlates  very  well  with  those  of 
Ishida's  two-dimensional  analysis.  The  worst  combination, 
i.e.,  the  combination  of  the  coarse  mesh  representation  and 
the  12-node  superelement  yields  the  F(a/b)  value  of  1.05 
which  is  about  10%  lower  than  Ishida's  solution. 

3.12.3  Rectangular  Prism  with  a Double  Edge  Crack 

Figure  10  shows  a rectangular  prism  of  width  2b,  thick- 
ness 2t  and  length  2L  containing  a double  edge  crack  of  depth 
a.  The  origin  of  the  global  coordinates  (X,Y,Z)  is  taken  at 
the  center  of  the  prism.  The  loading  condition  is  the  uniform- 
ly distributed  tension  load  applied  on  both  the  upper  and 
lower  surfaces  of  the  prism.  Two-dimensional  plane  strain 
solution  of  the  stress  intensity  factor  for  a double  edge 
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crack  in  a rectangular  plate  has  been  obtained  by  Bowie  and 
is  found  in  Reference  33. 

Since  the  geometry  and  the  loading  are  symmetric  with 
respect  to  the  plane  X=0,  Y=0  and  Z=0,  it  is  necessary  to  con- 
sider only  one  eighth  of  the  prism  for  finite  element  analysis. 
The  mesh  representations  used  in  this  analysis  are  the  same  as 
those  in  the  preceding  two  examples.  Symmetric  boundary  condi- 
tions are  imposed  on  the  plane  X=0,  Y=0  and  Z=0.  The  resulting 

1/2 

solution  are  shown  in  the  normalized  form  of  Y = K^/aa 

in  Table  17  and  in  Figure  14.  Two-dimensional  solutions  of 

1/2 

Bowie  are  expressed  as  Kj  = a(ira)  F(a/b),  where  a is  the 
applied  stress,  a is  the  crack  length  and  b is  the  half  plate 
width.  For  the  geometry  of  a/b  = 0.5,  L/b  -*■  °°,  his  computed 
value  of  F (a/b)  is  1.15,  where  L is  the  half  length  of  the 
plate.  When  the  fine  mesh  representation  is  used,  the  12-node 
superelement  adjacent  to  the  center  plane  gives  the  F(a/b) 
value  of  1.20  which  is  4.3%  higher  than  Bowie's  solution. 

When  the  coarse  mesh  representation  is  used,  the  resulting 
F(a/b)  value  in  the  20-node  superelement  adjacent  to  the 
center  plane  is  1.18  which  is  2.6%  higher  than  that  of  Bowie's 
analysis.  The  value  of  F(a/b)  obtained  by  the  combination  of 
the  12-node  superelement  and  the  coarse  mesh  representation  is 
1.06  and  is  7.8%  lower  than  Bowie's  value. 

The  above-mentioned  three  examples,  i.e.,  single  edge 
crack,  center  crack,  double  edge  crack  specimens  have  been 


analyzed  not  to  pursue  the  accurate  solutions  for  individual 
problems  but  only  to  confirm  the  possibility  of  the  three- 
dimensional  crack  element  based  on  the  assumed  stress  hybrid 
stress  finite  element  model.  The  35  and  85  elements  mesh 
representations  used  in  these  analyses  are  both  extremely 
coarse.  However,  the  resulting  stress  intensity  factor  solu- 
tions show  the  sufficient  engineering  accuracy. 

3.12.4  Compact  Tension  Specimen 

Compact  tension  specimen,  as  the  name  implies,  can  be 

made  more  compact  than  any  other  specimens  that  could  be  used 

for  fracture  toughness  KI  testing.  Due  to  its  compactness 

c 

the  specimen  has  become  widely  used  especially  in  the  fields 
where  economy  or  space  for  exposure  of  material  is  of  particu- 
lar interest.  Therefore,  the  analysis  of  compact  test  speci- 
men is  significant  to  the  application  and  development  of 
fracture  mechanics. 

Figure  15  shows  the  standard  compact  tension  specimen 
(ASTM  E-399-72) . The  geometry  is  a plate  of  thickness  t, 
width  W and  height  h with  a through  the  thickness  single  edge 
crack  of  uniform  depth  a.  Uniformly  distributed  line  load  of 
magnitude  T (force/length)  which  is  normal  to  the  crack  sur- 
faces is  applied  at  the  pin  holes  throughout  the  thickness. 
Test  results  obtained  by  using  compact  tension  specimen  are 
usually  analyzed  with  the  aid  of  the  two-dimensional  plane 
strain  solution  for  the  stress  intensity  factor.  However, 
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since  the  thickness  of  this  specimen  is  not  negligibly  small 
and  the  interior  points  are  not  remote  from  the  outer  surfaces, 
the  stress  field  in  the  interior  is  highly  three-dimensional 
and  hence  the  magnitude  of  the  stress  intensity  factor  appears 
to  vary  considerably  along  the  crack  front.  Therefore,  three- 
dimensional  analysis  is  essential  for  the  detailed  estimation 
of  the  specimen  behavior.  Several  authors,  Cruse  [46],  Tracey 
[59],  Barsoum  [67]  and  Yamamoto  and  Sumi  [70]  have  carried  out 
three-dimensional  analyses  of  this  specimen.  Because  of  the 
complexity  of  three-dimensional  crack  analysis,  numerical 
methods  have  been  used  by  such  authors.  Cruse  used  singular 
integral  method  (direct  potential  method) , Tracey  and  Barsoum 
performed  finite  element  analyses  based  on  assumed  displace- 
ment model.  Yamamoto  and  Sumi  applied  the  method  of  super- 
position of  analytical  and  finite  element  solution  to  deter- 
mine the  detailed  stress  field  of  this  specimen. 

The  orientation  of  the  axes  of  global  coordinates  (X,Y,Z) 
are  shown  in  Figure  15  with  the  origin  at  the  center  of  the 
crack  front.  Since  the  geometry  of  the  specimen  and  loading 
are  symmetric  across  the  plane  Y=0  and  Z=0,  only  one  quarter 
of  the  specimen  is  needed  in  the  finite  element  analysis.  In 
the  present  study,  the  pin  holes  are  removed  and  the  loading 
is  represented  by  a distributed  line  load  applied  along  the 
contact  lines  of  pins  and  pin  holes.  Twelve- node  upper  half 
superelement  which  consists  of  two  8-node  crack  front  elements 
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is  used  for  near-field  element,  while  8-node  hexahedron 
element  based  on  hybrid  stress  model  occupies  the  rest  of  the 
domain.  Two  different  specimen  geometries,  a/W  = 3/8  and 
a/W  = 5/8  are  considered.  In  both  cases,  the  Poisson's  ratio 
v = 0.3  is  used.  The  finite  element  mesh  subdivisions  for  two 
different  specimen  geometries  are  drawn  in  Figure  16  and  both 
models  have  85  elements,  168  nodal  points  and  504  degrees  of 
freedom.  In  these  finite  element  representations,  a quadrant 
of  the  specimen  is  divided  into  five  layers  at  equal  intervals 
by  planes  normal  to  the  crack  front.  Each  layer  contains  one 
superelement. 

For  the  specimen  of  a/W  = 3/8,  the  depth  of  the  crack 
contained  in  each  superelement  is  one  third  of  the  total  crack 
depth  a.  Along  the  X-axis,  the  nodal  points  are  arranged  at 
X/W  = -0.625,  -0.375,  -0.125,  0,  0.125,  0.375,  0.625 
Along  the  Z-axis,  the  nodes  are  located  at 
Z/W  = 0.125,  0.4,  0.6 

For  the  specimen  of  a/W  = 5/8,  the  depth  of  the  crack  in  each 
superelement  is  one  fifth  of  the  total  crack  depth  a.  The 
nodes  are  arranged  at  the  following  locations  along  the  X-axis 
X/W  = -0.875,  -0.625,  -0.375,  -0.125,  0,  0.125,  0.375 

I 

Along  the  Z-axis 
Z/W  = 0.125,  0.4,  0.6 

Symmetrical  boundary  conditions  are  imposed  on  the  plane  Y=0 
and  Z=0. 
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Figures  17  and  18  show  the  variation  of  the  magnitude  of 
the  stress  intensity  factor  thus  obtained  for  the  specimen 
of  a/W  = 3/8  and  a/W  = 5/8,  respectively.  They  are  also  listed 
in  Table  18.  The  corresponding  two  dimensional  plane  strain 
solutions  are  drawn  in  the  above  figures  for  convenience. 

These  solutions  are  normalized  in  the  form  of 


Y 


k /T ( 2W+a) 
Ixt  (W-a) 


where  T is  the  total  applied  load. 

The  resulting  solution  takes  the  maximum  value  in  the 
superelement  closest  to  the  center  of  the  specimen  and  gradu- 
ally decreases  as  the  free  surface  is  approached.  The  solu- 
tions obtained  by  Yamamoto  and  Sumi  are  plotted  in  Figures  17 
and  18  for  the  purpose  of  comparison.  For  the  geometry  of 
a/W  = 3/8,  the  correlation  between  these  two  solutions  is 
indeed  excellent.  For  a/W  = 5/S,  the  correlation  is  less  good 
than  that  of  a/W  = 3/8.  In  this  case  the  present  solution 
in  the  superelement  closest  to  the  center  of  the  specimen  is 
only  1%  above  the  plane  strain  solution,  while  Yamamoto's 
solution  is  5%  above  the  plane  strain  solution  at  the  center 
of  the  specimen.  It  is  seen  that  for  both  solutions  the  value 
of  K at  the  center  of  the  specimen  of  a/W  = 3/8  deviates  more 
• from  the  plane  strain  solution  than  that  of  the  specimen  of 

a/W  = 5/8. 

i '• 
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It  should  be  noted  that  for  the  present  solution  only  one 
set  of  matrix  equations  of  504  unknowns  is  solved  while  in  the 


method  by  Yamamoto  and  Sumi,  computation  of  finite  element 
solutions  should  be  performed  (2m  + 1)  times  where  m is  the 
number  of  finite  element  layers  across  the  thickness  of  the 
plate.  For  each  solution  they  considered  a mesh  pattern  with 
2106  degrees  of  freedom.  As  indicated  in  Table  19,  by  using 
a mesh  subdivision  with  fewer  degrees  of  freedom  than  that 
by  Tracey  and  Barsoum,  the  present  hybrid  stress  superelements 
provide  results  which  are  closer  to  the  solution  by  Yamamoto 
and  Sumi.  The  latter  are  used  here  as  reference  because  they 
were  obtained  by  finer  meshes  and  also  contain  a complete 
treatment  of  the  singular  behavior  at  the  crack  front. 

3.12.5  Cylindrical  Rod  with  an  Axisymmetrically  Located 
Penny  Shape  Crack  under  Uniform  Tension  Loading 

A cylinder  of  radius  R,  length  2L  containing  a penny- 
shaped crack  of  radius  a with  the  axis  of  symmetry  being 
normal  to  the  crack  plane  is  shown  in  Figure  19.  Cylindrical 
coordinates  (r,9,Z)  are  also  drawn  in  the  same  figure.  With 
reference  to  this  coordinates  system,  the  crack  occupies  the 
region  Z=0±,  0 < r £ a.  The  Z direction  is  the  longitudinal 
axis  of  the  cylinder  through  the  center  of  the  crack.  Uniform- 
ly distributed  tension  load  is  applied  on  both  upper  and  lower 
end  surfaces  of  the  cylinder.  The  semi-analytical  solution 
for  this  problem  has  been  obtained  by  Benthem  and  can  be  found 


in  Reference  36. 
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Since  the  geometry  and  loading  are  not  only  axisymmetric 
with  respect  to  the  Z-axis,  bat  also  symmetric  across  the  plane 
Z=0,  only  one  narrow  sector  of  the  cylinder  in  the  region  Z _>  0 
needs  to  be  modelled  for  finite  element  analysis.  Analysis  has 
been  carried  out  for  two  different  included  angles  9 = 10°  and 
9 = 5°,  and  for  two  different  finite  element  mesh  representa- 
tions, 396  degrees  of  freedom  and  810  degrees  of  freedom.  The 
396  d.o.f.  finite  element  model  is  shown  in  Figure  20(a).  It 
has  132  nodal  points  and  is  constructed  from  a twelve-node 
upper  half  superelement,  5 six-node  pentahedron  elements  and 
46  eight-node  hexahedron  elements  or  a total  of  52  finite  ele- 
ments. The  nodal  points  are  arranged  at  the  following  loca- 
tions along  the  r-axis 

r/R  = 0,  0.05,  0.1,  0.125,  0.15,  0.175,  0.2,  0.225,  0.25, 

0.275,  0.3,  0.35,  0.55,  1.0 

Along  the  Z-axis,  they  are  located  at 

Z/L  = 0,  0.05,  0.1,  0.2,  0.6,  1.0 

The  810  d.o.f.  finite  element  model  is  drawn  in  Figure  20(b). 
This  model  has  270  nodal  points  and  a twelve- node  upper  half 
superelement,  8 six-node  pentahedron  elements  and  107  eight- 
node  hexahedron  elements  or  a total  of  116  finite  elements. 

The  nodes  are  arranged  at  the  following  locations  along  the 
r-axis 

r/R  = 0,  0.05,  0.1,  0.125,  0.15,  0.175,  0.2,  0.225,  0.25, 

0.275,  0.3,  0.4,  0.7,  1.0 
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On  the  Z-axis,  they  are  at 

Z/L  = 0,  0.025,  0.05,  0.075,  0.1,  0.2,  0.4,  0.7,  1.0 
In  both  finite  element  representations,  one-eighth  of  the 
total  crack  radius  a is  included  in  the  superelement.  The 
far-field  elements,  i.e.,  6-node  pentahedron  and  8-node  hexa- 
hedron regular  elements  are  based  on  assumed  displacement 
finite  element  model  and  are  generated  by  using  isoparametric 
transformation  techniques.  To  simulate  the  axisymmetric 
behavior  of  the  cylinder,  all  nodes  are  fixed  in  the  circum- 
ferential direction.  The  nodes  which  are  ahead  of  the  crack 
plane  are  also  restrained  from  moving  in  the  Z-direction. 

The  stress  intensity  factor  K^.  thus  obtained  are  listed 
in  Table  20.  For  both  included  angles  the  fine  (116  elements) 
mesh  subdivision  yields  the  solution  within  1%  of  Benthem's 
result.  The  error  of  the  solution  obtained  by  using  the 
coarse  (52  elements)  mesh  representation  is  1.2%  for  0 = 10° 
and  2.9%  for  0=5°. 

In  order  to  check  the  influence  of  the  element  size  upon 
the  stress  intensity  factor  solution,  computation  has  been 
carried  out  for  the  finite  element  model  of  the  included  angle 
0 = 5°  by  using  the  crack  front  element  with  profile  edges 
two  times  larger  than  those  of  the  aforementioned  crack  front 
element.  The  fine  mesh  representation  (810  d.o.f.)  is  used. 
The  resulting  solution  is  = 1.617.  It  seems  that  the  effect 
of  the  element  size  is  not  severe  if  it  is  between  one-third 
and  one- tenth  of  the  crack  size  a. 
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Figure  21  shows  a cube  of  dimension  2b  x 2b  x 2b  contain- 
ing a penny-shaped  crack  of  radius  a(=  b/5) . The  global 
Cartesian  coordinates  (X,Y,Z)  and  cylindrical  coordinates 
(r,9,Z)  are  drawn  in  the  same  figure.  Uniformly  distributed 
tension  load  is  applied  on  both  top  and  bottom  surfaces  which 
are  parallel  to  the  crack  plane.  For  a penny-shaped  crack 
of  radius  a in  an  infinite  solid  under  a uniform  normal 
applied  stress  a at  infinity,  Sneddon  has  obtained  an  exact 
stress  intensity  factor  solution  KI  = 2o/a/ir  [32]  . 

Since  the  geometry  and  loading  is  symmetric  across  the 
planes  X=0,  Y=0  and  Z=0,  it  suffices  to  model  only  one-eighth 
of  the  cube  for  finite  element  analysis.  The  finite  element 
mesh  subdivision,  which  is  shown  in  Figure  24,  consists  of 
seven  planes  of  nodes  at  15°  intervals  about  the  Z-axis.  On 
these  planes,  the  nodes  are  arranged  at  the  following  dis- 
tances from  the  Z-axis. 

r/b  = 0,  0.05,  0.1,  0.15,  0.175,  0.2,  0.225,  0.25,  0.3,  0.5, 

1.0 

Along  the  line  drawn  from  the  crack  front  parallel  to  the 
Z-axis,  the  nodes  are  arranged  at 
Z/b  = 0,  0.025,  0.05,  0.1,  0.4,  1.0 

The  above  finite  element  representation  contains  284  nodes, 

852  degrees  of  freedom,  6 twelve-node  upper  half  superelements 
78  six-node  pentahedron  elements  and  138  eight-node  hexahedron 


elements  or  a total  of  222  finite  elements.  The  depth  of  the 
crack  included  in  each  superelement  is  one-eighth  of  the  crack 
radius.  All  far-field  elements,  i.e.,  6-node  pentahedron  and 
8-node  hexahedron  elements  are  based  on  the  assumed  displace- 
ment finite  element  model.  Symmetric  boundary  conditions  are 
imposed  on  the  plane  X=0,  Y=0  and  Z=0. 

The  resulting  stress  intensity  factor  solution  is 

Kj  = 1.58  or  Y * 0.990  for  all  crack  front  elements,  where 
1/2 

Y = KI/2a(a/ir)  ' . Since  the  Sneddon's  solution  for  a penny- 

shaped crack  of  radius  a in  an  infinite  solid  is 
1/2 

Kj  = 2a(a/ir)  ' , the  above  result  shows  - 1%  deviation  from  the 

Sneddon's  solution.  Tracey's  solution  for  a penny- shaped 


crack  in  a cylindrical  rod  is  shown  in  Figure  25  for  compari- 
son purposes. 

3.12.7  Semi-Circular  Surface  Crack  in  a Rectangular 


Prism 


Surface  crack  problems  occur  frequently  in  service  condi- 
tions. Figure  22  shows  a semi-circular  surface  crack  of 
radius  a which  lies  normal  to  the  surface  of  a rectangular 
prism  of  height  2b,  width  2b  and  thickness  b.  Radius  a of 
the  crack  is  one-fifth  of  the  half  width  b of  the  prism. 


The  global  Cartesian  coordinates  (X,Y,Z)  and  cylindrical 
coordinates  (r,9,Z)  are  also  drawn  in  the  above  figure.  The 


problem  is  to  obtain  the  variation  of  the  stress  intensity 
factor  along  the  semi-circular  crack  front  under  uniform 
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tension  loading  applied  at  both  upper  and  lower  surfaces  of 


the  prism. 

Since  the  geometry  and  loading  is  symmetric  across  the 
plane  X=0  and  Z=0,  only  a quadrant  of  the  prism  is  needed  for 
finite  element  analysis.  The  finite  element  mesh  representa- 
tion shown  in  Figure  24,  which  has  been  used  for  the  preced- 
ing example,  is  used  again  in  the  present  example.  Symmetric 
boundary  conditions  are  imposed  on  the  plane  X=0  and  Z=0. 

Figure  26  presents  the  variation  of  the  stress  intensity 
factor  thus  obtained  along  the  crack  front.  They  are 
normalized  by  Sneddon's  solution.  It  is  seen  that  Kj  takes 
the  minimum  value  in  the  superelement  of  the  deepest  penetra- 
tion while  the  maximum  stress  intensity  factor  occurs  in  the 
element  adjacent  to  the  free  surface.  Smith  [42]  obtained  the 
stress  intensity  factor  for  a semi-circular  surface  flaw  in 
a semi-infinite  solid  subjected  to  uniaxial  tension  by  using 
the  alternating  technique.  Tracey  solved  the  semi-circular 
surface  crack  in  a rod  of  semi-circular  cross  section.  Their 
results  are  also  shown  in  Figure  26  for  comparison  purposes.. 
The  present  solution  correlates  well  with  their  results. 

3.12.8  Quarter-Circular  Corner  Crack  in  a Rectangular 


Prism  m Tension 

Figure  23  shows  a quarter  circular  corner  crack  of 
radius  a in  a rectangular  prism  of  height  2b,  width  b and 


thickness  b.  Crack  radius  a is  one-fifth  of  the  width  of  the 


L. 


prism.  Uniformly  distributed  tension  loading  is  applied  on 
both  top  and  bottom  surfaces  of  the  prism.  The  global  coor- 
dinates (X, Y, Z)  and  (r,9,Z)  are  drawn  in  Figure  23. 

Since  the  geometry  and  loading  is  symmetric  across  the 
plane  Z=0,  one- half  of  the  prism  is  modelled  for  finite 
element  analysis.  The  same  mesh  subdivision  as  that  in  the 
preceding  two  examples  is  used  in  the  present  problem. 

The  variation  of  the  normalized  stress  intensity  factor 
solution  along  the  crack  front  is  presented  in  Figure  27. 
Also  shown  for  comparison  purposes  is  the  finite  element 
solution  obtained  by  Tracey  for  a quarter-circular  corner 
crack  in  a rod  of  quarter-circular  cross  section. 
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3.13  Concluding  Remarks 

A special  crack  front  element  has  been  developed  based 
on  the  hybrid  stress  finite  element  model  for  the  analyses  of 
arbitrarily  shaped  three-dimensional  crack  problems.  The 
superelement  is  compatible  with  most  of  the  existing  general 
purpose  finite  element  computer  programs.  The  CPU  time  for 
generating  a 12-node  crack  front  superelement  is  5.4  seconds 
on  an  IBM  S-370/165. 

For  obtaining  stress  intensity  factor  solutions  it  is 
not  necessary  to  use  any  interpolation  techniques  such  as 
C.O.D.  method,  virtual  crack  extension  method,  etc.,  which 
are  troublesome  particularly  for  a three-dimensional  crack 
problems  due  to  the  variability  of  stress  intensity  factors 
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along  the  crack  front.  In  the  present  crack  front  elements 
three  stress  intensity  factors  Kj  and  KIIX  are  unknown 

stress  parameters  and  hence  can  be  calculated  directly. 

The  displacement  compatibility  with  adjoing  far-field 
elements  is  completely  satisfied  across  interelement  bounda- 
ries whether  they  are  based  on  an  assumed  displacement  finite 
element  model  or  a hybrid  model.  However,  if  the  crack  front 
is  curved,  the  crack  front  elements  should  be  arranged  so 
that  the  crack  front  is  symmetrical  across  the  profile  face 
which  is  common  to  two  neighboring  crack  elements  in  order 
to  maintain  the  displacement  compatibility  on  this  face. 

The  stress  intensity  factors  are  assumed  to  be  constant 
within  each  crack  front  element.  To  determine  the  distribu- 
tions of  Kj,  and  along  the  crack  front,  a number  of 

special  elements  should  be  used  along  the  crack  front. 

In  a special  element  model  based  on  the  assumed  displace- 
ment method  a number  of  elements  must  be  used  around  each 
segment  of  the  crack  front  to  simulate  the  rapid  variation  of 
local  fields  with  9.  However,  in  the  present  method,  only  a 
single  superelement  is  needed  around  each  segment  of  the 
crack  front  and  thus  the  crack  front  region  can  be  modelled 

I 

with  much  fewer  degrees  of  freedom  compared  to  other  displace- 
ment-based crack  elements.  In  the  analysis  of  fracture  test 
specimens,  the  finite  element  mesh  representations  used  in 
the  present  study  contain  only  270  ~ 504  degrees  of  freedom 
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which  are  less  than  those  required  in  two-dimensional 
analyses  by  using  the  displacement-based  special  elements. 

For  three-dimensional  crack  problems  available  analyti- 
cal solution  is  limited  only  to  that  for  an  embedded  crack 
in  an  infinite  solid.  Due  to  this  lack  of  the  analytical 
solution  of  finite  domain  which  can  be  served  as  a basis  for 
comparison,  accuracy  of  the  present  method  is  still  an  open 
question  in  the  true  sense  of  the  word.  However,  the  solu- 
tions of  the  stress  intensity  factors  for  the  present  example 
problems  are  found  to  be  in  good  agreement  with  those  obtained 
by  other  semi-analytical  methods  and/or  two-dimensional 
analyses. 

As  the  present  crack  front  element  incorporates  a 
inverse  square  root  r stress  singularity  at  the  crack  front 
a priori,  it  appears  to  provide  no  additional  insight  into 
the  nature  of  singularity  at  the  intersection  point  between 
the  crack  front  and  the  free  surface.  Unfortunately,  at  the 
present  state  of  knowledge,  the  order  of  the  stress  singular- 
ity at  this  intersection  point  has  not  yet  been  known  and  no 

or 

numerical  technique  can  be  expected  to  provide  any  meaningful 
solution  for  this  fundamental  problem.  Therefore,  in  the 
present  study,  no  attempt  has  been  done  to  deal  with  the 
details  of  the  variation  of  the  stress  intensity  factor  solu- 
tion in  the  thin  boundary  layer  near  the  surface.  Even  for 
the  variation  of  the  stress  intensity  factor  approaching 


a free  surface,  there  exist  quite  contradictory  solutions. 

In  one  solution,  the  magnitude  of  increases  as  the  free 
surface  is  approached,  in  another  it  increases  and  rapidly 
decreases,  and  in  the  third  it  decreases.  However,  except 
in  a thin  boundary  layer  near  a free  surface,  the  present 
crack  front  element  appears  to  provide  highly  accurate  stress 
intensity  factor  solutions  all  along  the  crack  front. 

For  body  force  problems,  the  inconsistent  approach  for 
obtaining  equivalent  nodal  forces  is  adequate.  This  approach 
eliminates  the  non-uniqueness  problem  of  the  finite  element 
solution  for  the  type  A crack  front  element. 

The  present  method  can  be  easily  applied  to  thermal 
stress  problems.  The  temperature  variation  within  each 
element  can  be  expressed  in  terms  of  its  nodal  value  by  using 
interpolation  technique.  In  this  case  the  use  of  a 16-node 
crack  front  element  is  desirable  because  it  contains  constant 
strain  terms  in  the  boundary  displacement  assumption.  Present 
8-node  and  12-node  crack  front  elements  do  not  have  complete 
constant  strain  terms  and  so  may  provide  erroneous  results  in 
thermal  gradients. 

The  present  method  can  $lso  be  easily  extended  to  the 
determination  of  stress  intensity  factors  for  a crack  in  an 
anisotropic  material. 

.1 
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SECTION  4 

HYBRID  CRACK  ELEMENT  FOR  ANALYZING 
THROUGH  FLAWS  IN  PLATES  IN  BENDING 

4.1  Introduction 

The  characteristics  of  the  local  stress  and  displacement 
field  at  the  base  of  a through-the-thickness  stationary  crack 
in  a plate  has  been  analytically  investigated  for  both  symmet- 
rical and  antisymmetrical  out  of  plane  loadings  by  several  re- 
searchers. Williams  [101]  first  analyzed  this  problem  by  us- 
ing the  technique  of  Fadle  eigenfunction  expansions  based  on 
the  fourth  order  thin  plate  theory  of  Poisson  and  Kirchhoff 
and  found  that  the  elastic  bending  stresses  near  the  tip  of  a 
semi- infinite  straight  line  crack  vary  as  the  inverse  square 
root  of  the  radial  distance  from  the  crack  tip.  His  results 
were  incomplete  in  that  the  magnitude  of  the  local  stresses 
were  left  undetermined.  Later,  Sih  and  Rice  [103]  cleared  the 
way  of  finding  the  coefficients  in  the  eigenfunction  expansion 
through  an  application  of  the  theory  of  complex  functions. 

Sih,  Paris  and  Erdogan  [102]  modified  William's  results  to  de- 
fine stress  intensity  factors  in  bending  in  a manner  consist- 
ent with  Irwin's  definitions  for  extensional  problems  in 
plates . 

However,  the  above  solution  based  on  the  classical  plate 
theory  satisfies  the  stress  free  condition  along  a crack  sur- 
face only  in  an  approximate  manner.  The  three  physically 
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natural  boundary  conditions  of  vanishing  bending  moment, 
cwisting  moment  and  transverse  shear  stress  at  the  edge  of 
the  crack  are  contracted  into  two  conditions.  To  overcome 
this  shortcoming,  Knowls  and  Wang  [106]  formulated  the  problem 
in  terms  of  Reissner's  six  order  plate  theory  in  which  all 
three  boundary  conditions  on  the  crack  surfaces  can  be  satis- 
fied. They  reduced  the  formulation  to  singular  integral 
equations  and  obtained  the  form  of  the  bending  stress  distri- 
bution in  the  vicinity  of  a crack  tip  of  a vanishingly  thin 
infinite  plate  subjected  to  constant  bending  moment  at  infin- 
ity. They  found  that  the  bending  stresses  of  primary  impor- 
tance possess  singularities  at  a crack  tip  of  the  same  order 
as  in  the  classical  theory,  but  the  angular  distribution  of 
these  stresses  around  the  crack  tip  has  slight  discrepancies. 

Furthermore,  the  stress  resultants  Q and  Q become  infinite 

x y 

-3/2 

like  r ' in  the  classical  theory,  while  they  remain  finite 
as  r + 0 in  the  Reissner  theory.  Their  work  was  later 
extended  by  Hartranft  and  Sih  [107]  to  study  the  thickness 
effect  on  the  stress  distribution  around  the  crack.  Recently, 
in  the  framework  of  Reissner's  theory  Wang  [108]  obtained  the 
asymptotic  stress  field  for  a crack  in  an  infinite  plate  sub- 
jected to  constant  twisting  moment  at  infinity  by  using  the 
Fourier  integral  transform  techniques  to  reduce  the  problem 
to  a system  of  dual  integral  equations.  Accoring  to  his 
analysis,  the  singularity  of  the  stress  resultants  is  of  the 
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order  of  r in  case  of  twisting,  while  it  is  of  the  order 

-3/2 

of  r in  the  classical  theory. 


However,  it  is  noted  that  the  solutions  obtained  by  the 

aforementioned  methods  are  based  on  some  approximations  to 

a fully  three-dimensional  theory,  viz,  the  theory  of  bending 

of  a plate,  whether  it  is  a fourth  order  Kirchhoff  theory  or 

a sixth  order  Reissner  theory.  The  above  methods  assume 

linear  distribution  of  the  in-plane  stresses  a , a and  t 

c x y xy 

through  the  thickness  of  the  plate  and  thus  cannot  account 
for  the  nonlinear  disturbances  of  these  stresses  near  the 
crack  edges  and  the  surfaces  of  the  plate  which  are  important 
for  cracks  in  very  thick  plate.  For  these  cases,  the  use  of 
fully  three-dimensional  theory  is  necessary.  In  the  analyt- 
ical study  to  examine  the  effect  of  plate  thickness  on  the 
stress  distribution  around  the  crack  in  a plate  subjected  to 
out  of  plane  bending,  Sih  [119]  obtained  the  qualitative 
features  of  the  exact  solution  by  using  a three-dimensional 
asymptotic  expansion  of  stresses  and  displacements. 

As  described  above,  numerous  analytical  and/or  numerical,^ 
investigations  have  been  worked  out  for  the  flexural  problems 
of  plate  containing  a crack.  However,  very  few  reports  for 
the  application  of  finite  element  method  to  bending  problems 
of  cracked  plate  have  been  published  to  date,  though  a great 
number  of  plate  bending  finite  element  models  have  been 
developed  over  the  past  ten  years. 
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Wilson  and  Thompson  [109]  determined  the  bending  stress 


intensity  factor  for  a through  thickness  crack  in  a plate 
subject  to  purely  cylindrical  bending  by  using  conventional 
assumed  displacement  bending  elements  of  triangular  shape. 
Estimate  of  was  obtained  by  substituting  the  resulting 
values  of  displacements  at  different  points  near  the  crack 
tip  into  the  known  classical  solution  along  with  the  respec- 
tive coordinates  of  the  reference  point  and  solving  for  K^. 

The  resulting  values  appear  to  be  smooth  curves  when 
plotted  against  the  distance  of  the  reference  point  from  the 
crack  tip.  But  they  break  down  as  the  reference  point 
approaches  to  the  crack  tip.  The  break  down  is  due  to  the 
inability  of  the  assumed  displacement  function  to  represent 
the  singular  behavior  at  the  crack  tip.  The  best  estimate 
of  was  obtained  by  extrapolating  back  to  the  crack  tip 
from  the  results  calculated  at  several  points  just  beyond  the 
break  down  point  near  the  crack  tip. 

Luk  formulated  the  hybrid  stress  crack  element  based  on 
Reissner-type  plate  theory  without  implementation.  He  adopted 
the  second  scheme  described  in  Subsection  3.1  so  that  in  the 
vicinity  of  the  crack  tip  the  dominant  part  of  asymptotic 
solution  for  stress  couples  *nd  stress  resultants  obtained 
by  Hartranft,  Sih  [107]  and  Wang  [108],  is  included  in  elements 
surrounding  the  crack  tip.  In  Luk's  formulation  the  crack 
tip  superelement  consists  of  four  4-node  rectangular  basic 
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singular  elements  which  encircle  the  crack  tip.  The  assumed 
stress  field  in  the  interior  of  the  special  element  includes 
self-equilibrating  stress  terms  of  regular  polynomials  as 
well  as  the  known  singular  stress  solutions.  The  stress 
assumption  satisfies  a priori  the  traction  free  boundary 
condition  over  the  crack  surface.  In  the  boundary  displace- 
ment assumption,  both  the  transverse  deflection  w and  the 

1/2 

total  rotation  vary  as  r along  the  edges  which  adjoin  the 
other  basic  singular  elements.  However,  it  is  easy  to  point 
out  the  inconsistency  in  this  displacement  assumption.  Since 
the  exact  asymptotic  solution  for  the  displacement  field  is 
not  readily  available,  none  of  the  surface  integrals  in 
element  matrices  can  be  converted  to  line  integrals.  In 
constructing  a single  basic  special  element,  16  surface 
integrals  and  5 line  integrals  or  a total  of  21  distinct 
integrals  must  be  evaluated  for  generating  element  matrices. 
Some  of  the  integrand  involve  either  inverse  r singular- 
ity or  inverse  square  root  r singularity.  For  the  evaluation 
of  such  integrals,  special  treatment  is  required  to  obtain 
accurate  solutions.  Rhee  and  Atluri  [111]  used  the  hybrid 
stress  model  and  corrected  this  inconsistency.  They  solved 
several  pure  bending  problems  of  cracked  plates.  Before 
everything,  one  has  to  make  clear  the  asymptotic  behavior 
of  the  displacement  field  in  the  vicinity  of  the  crack  tip 
in  order  to  create  a crack  tip  element  model  for  practical 


use  based  on  the  2nd  scheme  of  the  hybrid  element  model  using 
Reissner's  6th  order  plate  theory. 

In  the  present  report,  a hybrid  crack  element  model  for 
bending  problems  is  developed  based  on  the  first  scheme,  using 
Poisson-Kirchhof f ' s 4th  order  plate  theory.  It  is  realized, 
however,  that  this  theory  is  only  applicable  to  thin  plates. 
Stress  and  displacement  fields  in  the  element  interior  are 
derived  by  the  complex  variable  method  and  conformal  mapping 
technique  and  they  satisfy  equilibrium  as  well  as  compatibil- 
ity conditions  in  the  general  boundary  value  problem  of 
cracked  plate.  The  procedure  to  obtain  the  superelement 
follows  closely  that  in  the  plane-element  given  in  Reference 
19.  In  later  chapters,  extension  to  bending  analysis  of  an 
isotropic  plate  with  wedge-shaped  notch,  an  anisotropic  plate 
with  through-the-thickness  crack,  and  bi-material  plate  with 
through-the-thickness  crack  located  parallel  or  normal  to  the 
interface  will  be  discussed. 

4.2  Stress  and  Displacement  Field  Equations 

The  asymptotic  behavior  of  the  Kirchhoff  stress-field  in 
the  vicinity  of  the  crack  tip  in  a plate  subject  to  out-of- 
plane bending,  which  has  been  obtained  by  Williams  [101]  and 
later  modified  by  Sih,  takes  the  following  form  when  expressed 
in  terms  of  stress  couples  in  Cartesian  coordinates  in  Figure 
28: 
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[-  Kx  cos|^  + K2  sin|^I  + 0 (r-1)  (r-KH) 


[Kx  sin-p-  + K2  cos^]  + 0 (r-1)  (4.1) 


where  and  K2  are  the  stress  intensity  factor  for  plate  bend- 
ing and  plate  shearing  respectively  (Figure  29) . The  two  coor- 
dinate system  (r,0,z),  (x,y,z)  are  defined  in  Figure  28,  h is 
the  plate  thickness  and  v is  Poisson's  ratio. 

The  corresponding  displacement  field  takes  the  following 

form: 
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(7+v)r3/2K1 
3/2  ( 3+v)  Gh 


(5+3v)K2 
3/2  ( 3 + V)  Gh 


tcos3i  . ao^I  coae, 


[sin§^  " 3wh\>)"  sir4] 


(4.2) 


where  G is  the  shear  modulus. 


4.3  Variational  Formulation 


The  following  notations  are  used  in  this  section. 


- stress  tensor 


- displacement  component 


- prescribed  displacement 


- component  of  body  force 


- surface  traction 


- prescribed  surface  traction  over  S 


- volume 


B (M  0 ) - complementary  energy  density 

Olp 

- stress  couples  per  unit  length 


- transverse  shear  stress  resultant 


- transverse  displacement  normal  to  the  midplane  of 
the  plate 

- deflection  along  3A 


- prescribed  deflection 


Hi 


- direction  cosine  of  the  normal  to  the  boundary 


v 

a 


n 


S 

s 


area  of  nth  element 

entire  boundary  of  nth  element 

portion  of  boundary  over  which  the  surface  traction 
is  prescribed 

portion  of  boundary  over  which  the  displacement  is 
prescribed 

interelement  boundary 
prescribed  stress  couple  over  CQ 

prescribed  transverse  shear  stress  resultant  over 
Ca 

boundary  of  V 

arc  length  along  boundary 


In  Hellinger-Reissner  variational  principle,  the  func- 
tional 77_  takes  the  form. 


I 

•< 


- E <"'v  B(<3ij>  + I °ij  <ui,j  + uj,i>  - FiuiIdv 

n n J j j ’ 


- ff  T.u.ds  - //  T. (u.  - u . ) ds} 
S 1 1 &u  1 1 1 

n n 


(4.3) 


The  independent  quantities  subject  to  variation  in  irR 

’ 

are  both  displacements  u^  and  stresses  . These  field 
variables  should  be  continuous  and  single-valued  in  each 
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element  and  compatible  along  interelement  boundaries. 
Specializing  above  functional  for  flexural  problem  of  plate 
and  excluding  body  force  for  simplicity  gives 


’R  - l lrUn  B(Mce)  - Mae“  c,6ldA 


+ (M  0v0w  - Q v w)ds 

c aB  B »a  a a 

an 


+ ;cu  IMc,sVW.c  ' ”,0.’  • Va(“  - “,lds} 


(4.4) 


where  a and  6 range  from  1 to  2. 

Introducing  the  displacement  compatibility  conditions 
along  the  interelement  boundaries  as  conditions  of  constraints 
through  the  Lagrange  multiplier  technique,  the  variational 
functional  is  modified  as 


%R  * l I'  B<*W  - Mo6u,aBldA 


+ (M  qvdw  - Q V w)ds 

ca  aS  B ,a  a a 

n 


+ /_  [M  v (w  - w ) - Q v (w  - w)  ] ds 
c^  ap  p ,ot  , ol  a a 

n 


+ (M  Rv  (w  - w ) - Q v (w  - w)  ] ds)  (4.5) 

c^  ap  p , a , a a a 


If  the  boundary  displacements  w and  w ^ are  so  chosen  that 
the  geometrical  boundary  conditions  are  satisfied. 


14 


w = w and  w = w along  S„  , the  functional  it  _ is 
» ot  / oi  un  mR 

rewritten  as 


%R  ' l - Ma6w,[.eldA 


+ V (Ma8v6”,a  - VaS)ds 

°n 


+ ;3An  ‘“aeW  - "a1  ' Vc(w  • *>Ids! 


(4.6) 


The  vanishing  of  <57rmR  with  respect  to  variations  of 
w and  w can  be  shown  to  give  the  following  Euler  equations 

(1)  Stress-strain  relation 

8B<*V  . _ 

w - ^ — — in  An 


(2)  Homogeneous  equilibrium  equation 


Ma6fOt6  0 


(3)  Compatibility  condition 


w = w,  w = w 

,a  , a 


(4)  Traction  boundary  condition 


M = M Q = Q 

a6  aB  a a 


(5)  Traction  equilibrium  condition 


M 0+  = M ",  Q + = Q ' 

aB  aB  a a 


in  An 


on  9An 


(4.7) 


on  C 


on  C. 
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' MtL'  i*"!  — 


where  the  superscripts  (+)  and  (-)  denote,  arbitrarily,  the 

left  and  right  sides  of  C as  C is  approached. 

n n 

4.4  Complex  Variable  Formulation 

It  is  recognized  that  the  complete  solution  of  plate 
bending  problem  can  be  obtained  from  the  equations  of  equilib- 
rium in  terms  of  displacements,  which  takes  the  form  of  a 
biharmonic  equation. 


fj 

■ 


V2V2w  = 0 


(4.8) 


The  moments  and  the  Kirchhoff’s  shear  forces  are  related  to 
the  displacement  in  the  following  manner. 


M 


M 


M 


xy 


-D(i%  + vii) 

9x * *y 

2 2 
9“w  . 9 w. 
-E)(v — j + — o 
9x^  dyZ 

9 2w 
9x9y 


= -D(l  - V) 


^ , 93w  , 93w 

-D  ( — + 5 


9x' 


-D<-^-  + 33” 


T 

9x  9y 


9x9y' 

i4 

9y' 


) 

) 


By  applying  the  method  of  complex  potential  originated  by 
Muskhelishvili  [112],  the  general  solution  to  the  Equation 
(4.8)  can  be  represented  in  terms  of  two  analytic  functions 
of  the  complex  variable  z(z  = x+iy) . 
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w(r,  9)  = Re  [ z<t>  ( z)  + x(z)l 


(4.9) 


where  z is  the  complex  conjugate  of  z. 

The  moments  and  Kirchhof f ' s shear  forces  corresponding  to  this 
displacement  field  are  expressed  in  terms  of  complex  potential 
<f>  (z)  and  x (z)  bY 

Mx  - -DRe  [2(1  + v)  4> ' (z)  + (1  - v)  {z<I>"  (z)  + x"  (z)  >1 


My  = -DRe  [2(1  + v)  4> ' (z)  - (1  - v)  {z$"  (z)  + x"  (z)  }] 


M = D ( 1 - v ) Im  [ z 4> " ( z ) + x"  ( z)  1 
xy 


Qx  = -4DRe[4>"  (z)  ] 


Qy  = 4DIm [$" (z) ] 


(4.10) 


or  in  combined  form 


M + M = -4D(1  + v)Re[4>'  (z)  ] 
x y 


My  - Mx  + 2iMxy  = 2D  (1  - v)  (z<f>"  (z)  +x"  (z)} 


Q - iQ  = -4D$"  (z) 
x y 


M«.ll) 


where  D is  the  flexural  rigidity  of  the  plate. 


D = Eh3/12 ( 1 - v2) 


v is  Poisson's  ratio.  ( )'  denotes  differentiation. 


The  derivatives  of  the  transverse  displacement  are  written  in 
the  form. 


3w 


= Re  [4>  (z)  + 2$'  (z)  + x'(z)] 


= Im  [4>  ( z)  ~ z^'  (z)  - x' (z) 1 


or  in  the  form 


(4.12) 


3w  . . 3w 


+ i-g^  = $(z)  + z4> ' (z)  + x'  (z) 


(4.13) 


where  ( ) denotes  the  complex  conjugate. 

The  stresses  and  displacements  given  in  Equations  (4.10), 

(4.9)  and  (4.12)  satisfy  exactly  the  compatibility  condition 
as  well  as  the  equation  of  equilibrium. 

4.5  Conformal  Mapping 

To  express  the  singularities  of  all  order,  the  following 
conformal  mapping  function  is  selected. 


z = to  U)  = 


(4.14) 


The  arguments  of  z and  £ are  limited  to  the  range  -tt<  arg  z 
<t\,  and  - arg  If/  for  a superelement  containing 

embedded  crack,  the  origin  of  the  coordinate  system  is  placed 
at  the  crack  tip  and  the  two  crack  surfaces  which  are  on  the 
negative  real  axis  of  the  z plane  are  mapped  to  the  imaginary 
axis  of  the  c plane  and  the  element  then  lies  on  the  region 
where  the  real  part  of  £ is  positive  [Figure  30] . 


Introducing  the  mapping  function  00(c)  into  Equations 
(4.10),  (4.9)  and  (4.12)  the  field  variables  are  given  by 

Mx  = -DRe  [2(1  + v)*'  (ci/w'  (C)  + (1  - v)  (ZJ7T T(*'  (C)/u' <?>>  ' 

+ (X' U)A>' (C)) ' }/u’ (5) 1 

My  = -DRe  [2(1  + v)*'  (O/w’  (O  - (1  - v)  {Z3TTT ( <I> * (?)/«'  (C)  ) ' 

+ (X*  (OA>'  U)  ) ' }/u'  (C)  ) 

Mvv  = D(l-  v)Im[{37cr(*'  (C)/m'(C))  ' 

xy 

+ (X'  (O/U'  (C) ) ' }/«*)'  (C)  ] 


Qx  = -4DRe[{$'  (C)Ao’  (?)  }'/u'  (C)  1 
Qy  = 4DIm[{$'  (C)A>'  (C)  }'/u'  (C)  ) 


(4.15) 


w = Re  [ 00(c)  4>  (C)  +xUH  (4.16) 

|^  = Re[$(c)  + {STTcJV  (5)  + X'  (0  )/«'  (C)  1 

|y  = lm[*(c)  - {Z3TcT«l» ' U)  +x'  (c)  }/u'  (c)  ] (4.17) 

The  Equations  (4.15)  and  (4.17)  may  be  combined  into  the  form 

M + M = -4D(1  + v)Re[*'  (C)/u' (C)  1 

x y 

M - M + 2iM  = 2D  ( 1 - v)  rZoTcT ( 4>  * (C)/w*  (C)  }' 

y x xy 
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+ {x' (?)/<*>'  (?)  }' 1/w'  (?) 


Q - iQv  = -4D{*'  U)/w' (?)  }'/<*>'  (?) 

x y 


+ iff  = *<?)  + u(?)f  (?)/o)''(?)'+  x'  (?)/<*>'  (?) 


(4.18) 


(4.19) 


The  complex  potential  <5  and  x which  are  analytic  functions  of 
z are  then  also  analytic  functions  of  ? and  can  always  be 
represented  by  polynomials  consisting  of  only  positive  powers 
of  5 There  are  two  alternative  approaches  to  the  derivation 
of  the  stress  and  displacement  fields  that  satisfy  the  stress 
free  boundary  condition  over  the  crack  surface. 

Approach  1 

In  the  first  approach,  4>(?)  and  x(?)  are  assumed  indepen- 
dently and  all  the  field  variables  M , M , M , Q , Q , W,  W 

x y xy  x y / x 

and  W are  expressed  in  terms  of  ?.  Then  they  are  substi- 
r y 

tuted  into  the  stress  free  boundary  condition  over  the  crack 
surface  to  obtain  the  relations  between  the  coefficients  of 
$(?)  and  x(?)-  The  potential  functions  $(?)  and  x(?)  are 
assumed  as 


2N  i 
*(?)  = £ b.?11 

j=i  1 


I (Bj  + iB2N+j) ? 


= j^1 1 (62j-l  + lS2N+2j-l) 5 ^ + ( B2j  + 102N+2j * 5 ^ 1 


2N+2  • 

X(C)  = E c.£3 

j-3  3 
2N+2 


= Z (Yj+iY2N+j)C 


j=3 

^ 2i+l 

= 2 I (Y^^xi  + ^ Y OXTJ.  0 4 J.1  )r3  ■L  + 

j 


2j+2 

^ 1 ' r2j+l  T xr2N+2j+l; ^ T (Y2j+2  + 1Y2N+2j+2) C 1 


(4.20) 


in  which  N is  a finite  integer  and  B's  and  y's  are  real 
constants. 

Taking  as  boundary  conditions  along  the  crack  surface 
(C  = in)»  the  usual  Kirchhoff  condition  for  a free  edge;  viz.. 


My  = 0 


(4.21) 


vy  s 0 


where  V..  = 


(4.22) 


3M 

Q + — ^ 

y 


3x 


Over  the  crack  surface.  My  and  Vy  are  expressed  as 


N 


My  = -D  E (-1)  3“1[[{(l  + v)  (2j  - 1)  - |(1  - v)  (2j  - 1)  (2j  - 3)}62N+2j 


- ±(1  - v)  ( 2 j + 1)  ( 2 j - 1)  y 


2N+2j+l 


1 n 


2 j~3 


+ [{(l  + v)2j  -i(l- v)2j(2j  - 2)}&2j 
- i(l-  V)  (2 j + 2)2jY2j  + 2]n2j"2] 
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(4.23) 


B*.-:--  -t 


_ N . , 

V,  = - £ (-1)-,”-l[[{2(5-  v)  (2j  - 1)  (2j  - 3) 

+ (1- v)  (2j  - 1)  (2j  - 3)  (2j  - 5)}B2j_1 
+ (1-  v)  ( 2 j + 1)  ( 2 j - 1)  ( 2j  - 3)Y2j+1]n2j"5 


- [{2(5  — v) 2j (2j  - 2)  + ( 1 — v) 2 j ( 2 j - 2)  ( 2 j - 4) }B2N+2j 
+ (1- v)  (2j  + 2)2j(2j  - 2)Y2N+2j  + 2]n2j'4]  (4.24) 


Substituting  Equation  (4.23)  and  (4.24)  into  Kirchhoff  stress 
free  condition  (4.21)  and  (4.22),  all  y's  can  be  expressed  in 
terms  of  B's. 

Y2j+1  * 

Y2j+2  = <J<23»e2j 

Y2N+2j+l  = q(22  " ^2N+2j-l 

Y2N+2j  + 2 = P<22>  S2N+2j 
where 

+1) 

q ( j ) - 8/  ( j + 2)  (1-  v)  - 1 

Hence,  all  the  field  variables  M , M , M , Q , Q , w,  w 

x y xy  x y / x 

and  w y can  be  expressed  in  terms  of  6's  only. 

The  moments  are 


(4.25) 
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I 


M = E a . 8 . 
x i=l  1 1 


My  '^i8! 

4N 

M = E c.g. 
xy  i=1  ipi 


where 


(4.26) 


a2j-1  = -D[{(1  + v)  +|(1- v)  (2j  + l)p(2j  - 1)  }(2j  - l)Re(c2j"3) 

+ - v)  ( 2 j - 1)  ( 2 j - 3)Re(c2i;2j_5)  ] 

a2j  - -D[{(1  + v)  +J(1-  v)  (2j  + 2)  q (2 j ) }ReU2j'2) 

+ J(l-  v)2j(2j  - 2)ReuVj_4)  ] 

a2N+2j-l  = °H(1  + v)  + J(1-  v)  ( 2j  + 1)  q ( 2 j - 1)  }(2j  - l)Im(C2j“3) 
+ J(l-  v)  (2j  - 1)  (2 j - 3)Im(?Vj“5)  ] 
a2N+2j  = D({<1  + v>  +Jd- v)  ( 2 j + 2)p(2j)  }2jImU2j~2) 

+ J(l-  v)2j(2j  - 2)Im(52C2j“4)] 
b2j_1  = -D[{(l  + v)  - J(l-  v)  ( 2 j + 1)  p (2j  - 1)  } (2 j -l)Re(c2j“3) 

- j(l- v)  (2j  - 1)  (2j  - 3)Re(?Vj-5)] 
b2j  - -D[{(1  + v)  - J(l-  v)  (2j  + 2)  q (2  j ) }2jReU2j_2) 

- \(1- v)2j(2j  - 2)ReU2£2j"4)  ] 

b2N+2  j-1  = Df{(1  + v)  - J(l-  v)  (2j  + l)q(2j  - 1)  }(2j  - l)Im(C2j"3) 


t* 


- \a  - v)  ( 2 j - 1)  ( 2 j - 3)lm(? Vj~5)  ] 
b2N+2 j = DIU1  + v)  - i(l-  v)  (2j  + 2)  p(2j  ) }2jlm(^2j"2) 

+ |(1-  v)2j(2j  - 2)Im(?Vj“4)  ] 

c2j-i  " Jd(I-v)  (2j-l)  [(2j-3)Im(?2c2j_5)  + (2j+l)p(2j-l)Im(C2j“3)J 
c2j  = (1-v)  2 j [ (2 j-2)  Im (T2?2^-4)  + (2j+2)q(2j)Im(£2^-2)  ] 

c2N+2j-l  =iD(1“v)(2j'1)[(2j“3)Re  (^2^2j"5)  + (2j+l)q(2j-l)Re(C2j“^] 
C2N+2 j =iD(1_  v)2j  [(2j-2)Re(?2C2j_4)  + (2j+2)p(2j)Re(C2j'2)  ] 

(j  = 1,  2,  3, , N-l,  N) 


The  Kirchhoff's  shear  stresses  are 


2N 

Qx  = j (dl j 0 j + d2j62N+j) 
2N 

Qy  *j=1Cd2j6j  ~ dljS2N+j) 


where 

dxj  = -Dj ( j - 2) Re ( ? j-4) 

d2j  = Dj(j  - 2)Im(?j_4) 
The  displacement  w is 
4N 


where 


(4.27) 


(4.28) 
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f2j-l  = Re^2^2j_1)  +P(2j  - l)Re(c2j+1) 
f2j  = ReU2C2j)  +q(2j)Re(C2j  + 2) 
f2N+2j-l  = -Im(?2^2j-1)-q(2j-l)Im(C2j+1) 


f2N+2j  = “ImU2;2j)  -p(2j)lm(c2j+2) 


(j  = 1,  2, 


The  slopes  w and  w are 

# x 9 y 


w X = Z gi8i 

'x  i=l  1 1 


w-y 


where 


(4.29) 


g2j-l  = 11  +i(2j  + DP(2j  - 1)  ]ReU2j-1)  +i(2j  - 1)  Re  (cVj~3) 
g2j  = d + (j  + l)q(2j)  JRe(;2j)  + jRe(?2C2j_2) 

g2N+2 j-1  = “ [1  + ( 2 j + 1)  q ( 2 j - 1)  ] Im  ( £ 2 j_1)  - ^ ( 2 j - 1)  Im (?V j_3) 

g 2N+  2 j = -I1+  <3  + 1)  P (2  j ) ]lm(?2j)  - jlm(^2?2j"2) 

h2j-1  = [1  - ^ (2j  + 1) p (2j  - 1) ]lm(c2j"1)  - i(2j  - 1)  Im  ( cV  j~3) 

fc2j  » [1-  (j  + l)q(2j)]Im(c2j)  - jIm(c2C2j"2) 


h2N+2j-l=  +Dq(2j  - DJReU2^1)  -|(2j  - 1)  Re  (cVj“3) 

h2N+2j  = [1  " (j  + l)p(2j)  ]Re(C2j)  - jReU2^2j"2) 

(j  = i,  2,-“,  N) 
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Approach  2 

In  the  second  approach,  first  only  $(;)  is  assumed  in 
terms  of  5.  The  stress  free  boundary  condition  is  expressed 
in  an  integral  form  by  using  the  Lechnitzky's  extension  of 
Muschelisvili ' s method  and  x'(£)  is  related  to  $(C) . 

The  stress  free  boundary  condition  over  the  crack  sur- 
face can  be  written  in  an  integral  form  as 


n<f>  ( z)  + z$'  (z)  + x'  (z)  = 0 


(4.30) 


where 


n = - (3  + v)/(l  - v) 


(4.31) 


For  the  following  derivation,  it  is  expedient  to  introduce  a 
new  variable  ^(z)  such  that 


*(z,  = 


(4.32) 


Using  the  mapping  z = uj(C)  = ^ and  the  new  variable  (z)  , 
Equation  (4.3g^  is  written  as 


n$  (O  + u (?)' 


+ 1J3  (Q)  * 0 


(4.33) 


Equation  (4.33)  is  satisfied  by  taking  the  complex  potential 
\p  (c ) in  the  form 


<MO  » -n*(-c)  -w(-c)*' (c)/w' (c) 


(4.34) 


The  potential  4>(C)  can  be  assumed  in  the  form  of  a polynomial 
C as 


N 

4>U)  = I b.  S3  (4.35) 

j*l  D 

where  bj  = + ^®N+j'  N is  a finite  integer  and  8's  are  real 

constants.  In  case  of  in-plane  elasticity,  SN+2  are  zero, 
because  the  corresponding  terms  of  4>(5)  and  iMC)  contribute 
no  contribution  to  stresses.  However,  in  flexural  problem, 

8n+2  is  not  necessarily  zero  and  cannot  be  eliminated  a priori. 
Inserting  Equation  (4.35)  into  Equation  (4.34)  gives 

N 

4>U)  - - E [nb.  (-I)3  + h>.  R3  (4.36) 

j=l  3 * J 


From  Equation  (4.18),  (4.16),  (4.19),  (4.35)  and  (4.36),  one 
can  write 


where 


Aj  * -Re(j;j"2) 


j . ,i- l! . ^ . iUsi-: 


n _ , ri  - 2 c2  A .(-l)3 

BN+j  ' 11  4 p + n_l 4l3C 


-j(j  ' 2)C 


j-4 


’N+j 


-ij  (j  - 2)  C 


j-4 


Ej  - Re[C2?jl  - y~-j[2n(-l)j  - j]RetCj  + 

fN+j  = -Imt?2Cj]  - y-i-j[2n(-l)  j + j]Im[C 

g.  * S?  - n(-l)  + ij^”2Im  [£2] 

J ur 

?N+j  = i{^  “ + j?j~2RetC2]} 

for  j * 1,  2,  3, , N - 1,  N 


It  is  obvious  that  these  two  approaches  will  yield  identical 
results.  The  latter  appears  to  be  expressed  in  simpler  form 
but  the  first  one  is  more  straightforward.  From  the  computa- 
tional point  of  view,  however,  the  first  approach  is  more 
efficient  and  convenient  than  the  second  approach  especially 
for  the  construction  of  the  superelement  for  flexural  problems. 

4.6  Finite  Element  Formulation 

The  stresses  and  displacements  given  in  Equations  (4.26), 
(4.27),  (4.28)  and  (4.29)  or  Equations  (4.37),  (4.38)  and 

(4.39)  satisfy  equations  of  equilibrium  and  compatibility  and 
kinematic  boundary  conditions  in  the  general  boundary  value 
problem  of  a plate  containing  a crack.  These  solutions  for 
the  stresses  and  displacements  are  used  to  construct  around 
the  crack  tip  a superelement  which  accounts  for  singularities 
of  all  order  under  Poisson-Kirchhof f theory.  For  elements 
far  away  from  the  crack  front,  the  regular  hybrid  and/or 
conventional  assumed  displacement  elements  can  be  used. 
Compatibility  of  displacements  and  equilibrium  of  tractions 
between  the  crack  tip  superelement  and  the  surrounding  regular 
elements  is  maintained  through  a Lagrange  multiplier  technique 
based  on  the  modified  Hellinger-Reissner  variational  principle 
[Eq.  (4.6)]. 


Since  some  of  the  Euler  equations,  i.e.,  w 


3B«V 


,a8 


DM 


a 8 


and  M 0 o=0  are  satisfied  a priori  in  the  crack  tip  super- 

ap,  dp 

element,  they  can  be  removed  from  the  modified  Hellinger- 
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Reissner  functional  [Eq.  (4.6)].  Thus,  the  functional 
corresponding  to  the  superelement  which  is  designated  as  Aj^ 
becomes 


it  = (-M  „v„w  + Q v w)dS 

mR  'aA^  a 6 8 ,a  a a 


- 4 . (-M  „v„w  + Q v w)dS 

2 aA^  a6  B ,a  a a 


+ fC  (Ma6V,a  “ QaV)dS 

al 


(4.40) 


When  the  assumed  stress  hybrid  elements  are  used  for  far-field 
elements,  the  homogeneous  equilibrium  equation  Mag>ag  = 0 is 
satisfied  a priori  in  these  elements.  Thus,  the  functional 
for  far-field  elements  is  written  as 


m 


mR  * £,(/a  -B(MaS,dA  + + W*>dS 

n-  fc  n 


n 


+ rc  (Mc.bV,c,  - Vc*,dS 

°n 


(4.41) 


Therefore,  the  total  hybrid  function  irmR  becomes 

TT  = TTC_  + TTRr 

mR  mR  mR 


(4.42) 


As  mentioned  before,  the  far-field  elements,  n - 2,  . . . m, 
can  be  simply  generated  through  the  conventional  assumed 
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displacement  finite  element  model  provided  the  displacement 
field  is  the  same  along  the  interface  between  these  elements 
and  the  crack  tip  super element . 

4.7  Derivation  of  Element  Stiffness  Matrix 

The  above  functional  tt  can  be  expressed  in  a matrix  form 
as 


* * ;3a,!T“  ds  ' 1 /sa1!ts  as  - rc  fa  ds 

1 1 

NE  -i  rp  m m 

+ Z [±  /_  o Co  dA  + / Tiu  ds  - / T n dS] 

n=2  1 An d n ~ ~ 

n 


where 


T = ■(  M 


Q v 

+ Q V 

X X 

y y 

M v 

+ M V 

X X 

xy  y 

M v 

+ M v 

xy  x 

y y 

-w  v + w v 
,v  x ,s  y 

-W  V - W V 

, v y , s x 


(4.43) 


This  functional  is  analogous  to  that  used  by  Tong,  et  al.  [19] 
who  treated  the  plane  stress  crack  problem  in  two  dimensional 
elasticity.  In  the  following  derivations  of  the  crack  tip 
superelement,  it  is  necessary  to  consider  only  the  first 
three  integrals  in  Equation  (4.43). 

Using  Equations  (4.26),  (4.27),  (4.28)  and  (4.29)  or 

Equations  (4.37),  (4.38)  and  (4.39)  the  assumed  stress  and 

displacement  modes  can  be  expressed  in  terms  of  undetermined 
real  parameter  B ' s as 

T = RB  (4.44) 

u = UB  (4.45) 


The  assumed  boundary  displacements  u can  be  expressed  in  terms 
of  the  generalized  nodal  displacements  as 


u = Lq 


(4.46) 


where  L is  the  one-dimensional  interpolation  function  matrix 
defined  on  the  element  boundaries  excluding  the  crack  surface. 
The  generalized  nodal  displacement  q can  be  selected  in  the 
form 


Si  * {w'  w,x'  w,y>i 


(4.47) 


The  interpolation  function  must  be  chosen  such  that  when  the 
nodal  displacements  for  the  neighboring  elements  are  matched, 
the  displacements  along  the  corresponding  interelement  bound- 
ary are  the  same.  We  can  use  the  osculatory  Hermite  poly- 
nomials (s)  to  represent  the  assumed  boundary  displace- 

ment w in  terms  of  generalized  nodal  displacement  q.  For 
instance  along  the  edge  AB  in  Figure  32,  w is  expressed  as 
follows 


w = H<»  <S)WA  + (S)WB  ♦ H<«  (slw^  + (s)»  Sb 

= H<»  (S|WA  ♦ H<»  (S)WB  + (s)  |-.y»  * V,y> 


* 


H10  (s)(-vw  +vw  ) 
2 y ,x  X , y 


(4.48) 


where 


Hm<«>  - 1 - 3(f>2  + 2(f)3 

H^’(s)  = 3(f)2  - 2(f)3 

(4.49) 

»u  (=)  - ‘if  - 2<f>2  + 'f)3) 

“S'  <•)  * -* ((f) 2 - <!>3i 


i is  the  distance  between  nodes  A and  B,  and  s is  measured 
from  A. 

From  Equation  (4.48),  we  obtain 


dH01}  dE02  dHU> 

w = — T^-w,  + — ■ ~W—  + — -j— — (— V w + V w ) 
,s  ds  A ds  B ds  y , xA  x ,yA 


dH 


(1) 


12 


a—  -(-V  W + v w ) 
ds  y ,xB  x ,yB 


(4.50 


The  normal  slope  along  the  edge  AB  must  be  linear  in  s to 
maintain  normal  slope  compatibility  along  the  entire  edge. 
Hence,  we  may  use  a linear  interpolation  function  to  repre- 
sent it  as 


(1 


(4.51 


The  assumed  boundary  displacement  u can  thus  be  expressed  in 
terms  of  the  generalized  nodal  displacement  q.  This  leads  to 


to  a cubic  variation  for  w and  a linear  variation  for  the 
normal  slope  on  the  boundary.  The  displacement  interpolation 
function  matrix  L is  given  in  Table  24.  It  is  noted  that  the 
boundary  displacement  u are  assumed  only  on  interelement 
boundaries  and  it  is  not  necessary  to  assume  those  over  the 
crack  surfaces. 

Substitution  of  Equations  (4.44),  (4.45)  and  (4.46)  into 

the  first  three  terms  of  Equation  (4.43)  yields  equation  of 
the  form 

7 r = 6TGq  - -|bTH6  (4.52) 

H are  given  by 

(4.53) 

? = I /3A  (?Ty  + HT?)ds  (4.54) 

The  stationary  condition  of  the  functional  it  with  respect  to 
variations  of  3 yields 

Gq  - HB  = 0 (4.55) 

By  solving  for  S from  Equation  (4.55)  and  substituting  back 
into  Equation  (4.52),  the  functional  n can  be  expressed  in 
terms  of  the  generalized  nodal  displacements  q only,  i.e. 


The  matrices  G and 


5 * 
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1 

TT  = 2? 


(4.56) 


where  K = GTH~1G  (4.57) 

K is  the  stiffness  matrix  of  the  crack  tip  superelement. 

It  is  well  known  that  for  a successful  implementation 
of  this  finite  element  model  the  number  of  B's  should  be 
greater  than  or  at  least  equal  to  the  number  of  degrees  of 
freedom  minus  the  number  of  rigid  body  modes  for  each  element. 
This  condition  is  understood  either  as  the  necessary  condition 
to  make  the  element  kinematically  stable  or  as  the  solvability 
condition  of  q's  in  terms  of  B's. 

4 . 8 Geometry  of  the  Crack  Tip  Superelement 

The  superelement  can  be  of  arbitrarily  polygonal  shape 
except  that  it  needs  to  contain  a crack  [Figure  31].  The 
vertices  of  the  polygon  are  chosen  as  the  nodal  points.  If 
the  outer  boundary  of  the  superelement  is  selected  as  a 
regular  polygon,  the  integration  of  element  matrices  G and  H 
is  simplified.  A further  numerical  simplification  can  be 
achieved  if  the  superelement  contains  a crack  along  an  axis 
of  symmetry  of  the  polygon,  since  certain  symmetrical  and 
antisymmetrical  properties  can  be  used  in  the  computation  of 
element  matrices.  It  should  be  pointed  out  that  the  nodal 
point  may  not  be  located  on  the  crack  tip.  since  any  terms 
which  involve  the  displacement  at  the  crack  tip  are  cancelled 


-i  — rzz 


in  the  integrations  of  the  element  matrices  G and  H which  are 
taken  all  around  the  element  boundary. 

For  the  present  formulation,  the  superelement  generated 
is  a nine-node  rectangular  shape  element  with  27  degrees  of 
freedom  as  shown  in  Figure  32.  In  this  case  at  least  24  B's 
must  be  used  in  the  representation  of  stress  and  displacement 
field  to  fulfill  the  condition  described  in  the  previous 
section. 

4.9  The  Integration  of  the  Matrix  G and  H 

It  is  seen  that  in  constructing  the  crack  tip  superele- 
ment only  two  line  integrals  along  the  element  boundary  need 
be  evaluated.  On  the  basis  of  the  Poisson-Kirchhof f theory 
of  thin  plates,  the  stress  and  displacement  field  given  in 
Equations  (4.44)  and  (4.45)  satisfies  the  stress-free  condi- 
tion along  the  crack  surface  only  in  an  approximate  manner. 
The  three  physically  natural  boundary  conditions  of  prescrib- 
ing bending  moment,  twisting  moment  and  transverse  shear 
stress  are  replaced  by  the  following  two  conditions: 

Mv  = 0 

on  crack  surface 

Q + -iir  - 0 

Therefore,  neither  the  shear  stress  Q nor  the  twisting  moment 

M^g  is  identically  zero  on  the  crack  surface.  The  line 
. t 

integral  of  T u over  the  crack  surface  can  thus  be  written 


1 I S ds  * WQu  + Mxy",x)dx  + ;A0(Qu  - Mxyw,x)dx 


3M  n 

= ‘W0  + -^)wdx  + V? 


9M  n a 

+ 'ao(0  + ^>"dx  - Mxyw|A  - Mxyw  1 P 


(4.58) 


in  which  0 denotes  the  crack  tip  and  A and  F are  the  nodes  at 

the  intersection  of  the  crack  surface  and  the  side  of  the 

crack  tip  super element.  It  should  be  noted  that  in  the  above 

expression  the  values  of  M are  those  on  the  crack  surface 

xy 

side  of  these  nodes. 

T 

In  a similar  manner,  the  line  integral  of  T u over  the 
crack  surface  becomes 


/„_.TTu  ds  = M ‘w|a 
FOA~  - xy  1 f 


(4.59) 


These  terms  are  represented  by  concentrated  loads  at  the 
nodes  and  should  be  added  to  the  matrices  G and  H.  Thus,  it 
is  seen  that  the  boundary  displacements  w and  w need  not 

9 v 

be  assumed  over  the  crack  surface.  The  process  for  integrat- 
ing the  G and  H matrices  does  not  involve  any  singular  terms, 
hence,  can  be  carried  out  numerically  without  any  special 
treatment. 

When  the  first  approach  is  adopted  to  obtain  stress  and 
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displacement  fields,  all  field  variables  M , M , M , Q , Q , 

x y xy  x y 


w,  w „ and  w „ are  expressed  explicitly  in  terms  of  B's. 
# x # y 


Therefore,  it  is  a straightforward  procedure  to  calculate  all 
the  terms  of  the  matrices  R and  U and  then  matrices  G and  H 
are  computed  according  to  Equations  (4.53)  and  (4.54)  by 
standard  numerical  integration  technique.  However,  if  the 
second  approach  is  adopted  the  field  variables  are  not 
expressed  distinctly  in  terms  of  B's  but  are  given  in  the 
form  of  complex  stress  and  complex  displacement.  In  this 
case,  the  computation  of  the  matrices  G and  H are  accelerated 
in  the  following  manner.  From  the  first  two  equations  of 
Equations  (4.37),  one  can  easily  obtain 


2N 


M*  "iMw  = D 1 P-iU.OB., -p(C,  s;  8^) 
A j =i  j j j 


(4.60 


2N 

Mv  + iMvv  = D E hi  U' S)  = hU'  S S 
y xy  j=1  3 3 


where 

Pj  = (1  + vjAj  - (1  - vjBj 

hj  = (1  + v)  A . + (1  - v)  B . 

T T~ 

One  can  then  express  T u and  Turn  the  form 


(4.61 


T 

T u = -Im((ipv  + hv  )g  - icvf] 


(4.62 


T u = -Im[ (ipv  + hv  ) (w  + iw  ) - icvtf]  (4.63) 

~ ~ x y,x  ,y 


where  v is  the  complex  direction  cosine  defined  by 


v = v + IV 
x y 


P,  h,  g,  c,  and  f are  complex  functions  of  x and  y defined 

in  Equations  (4.60),  (4.61),  (4.39),  (4.37)  and  (4.38).  w, 

w , w are  assumed  boundary  displacements  given  in  Table  24, 
t x t y 

The  integration  for  the  matrices  H and  G can  now  be 
can  now  be  carried  out  by  standard  numerical  schemes.  For 
instance,  the  integrations  of  H and  G for  the  particular 
geometry  of  Figure  31  become 


6TH6  = /3A  TTu  ds 


= "Re^AB  + /CD  + /EFVx(pg"Cf)dS] 


- + /DEvy (hg  + cf)ds]  + Im[hf ] 


(4.64) 


8TGq  = /3A  TTu  ds 


"RS[/AB  + fCD  + ;EFVx{p(w,x  + iw,y>  ’ C^}dS] 


Im[fBC  + -fDEvy(h(wfX  + iw  y)  + cw}ds]  + Im[hw]  (4.65) 


in  which  pg,  hg,  cf  and  hf  are  expressed  in  the  form 


pg 


2N  2N 

DJ1J16»6n,P.,ntV.)/2 

m=ln=l 


hg 


2N  2N 

DIE  BJ  (hg  + h g )/2 
m n m3n  n3m  ' 

m=ln=l 


2N  2N 
cf  = D Z Z 
m=ln=l 


BJJc  f + cf)  / 2 
m n m n n m 


hf  = 


2N 
D Z 


2N 

Z 


m=ln=l 


Bft  (h  f + hfJ/2 
m n m n n m 


(4.66) 


When  the  crack  lies  along  the  x-axis,  the  first  2N  (for 
approach  1)  or  N (for  approach  2)  terms  correspond  to  a 
symmetrical  solution  and  the  last  2N  or  N terms  correspond 
to  antisymmetrical  solution  with  respect  to  the  x-axis. 
Therefore,  the  matrix  H has  the  property 


Hj,M+j  HM+j,j  0 


(4.67) 


where  M = 2N  for  approach  1 
M = N for  approach  2 

Thus,  H takes  the  form 


H 

2M*2M 


?1 

MxM 


?2 

MxM 


(4.68) 


This  property  reduces  the  effort  for  the  computation  of  H 
matrix  and  its  inverse. 
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In  plane  stress  crack  problems,  the  H matrix  is  reduced 
to  (M-l)  x (M-l)  matrix  by  eliminating  BM+2,  but  in  flexural 
problems,  the  size  of  the  H matrix  can  not  be  reduced. 

For  the  crack  tip  superelement  in  Figure  3 2,  the  stiff- 
ness matrix  is  given  in  Table  29.  In  this  example,  the  mate- 
rial constants  E=l,  v=0.3  are  used.  For  the  assumption  of 
field  variables  24  B's  are  used  and  five-point  Gaussian  quadra- 
ture is  chosen  for  the  integration  of  the  G and  H matrices. 

4.10  Stress  Intensity  Factor 

The  bending  stress  intensity  factors  and  K2  on  the 
basis  of  the  Kirchhoff's  hypothesis  are  originally  defined  by 
Williams  [101].  But  later  Sih  et  al.  redefined  them  consider- 
ing the  similarity  between  out  of  plane  bending  and  in  plane 
stretching  problems.  However,  the  physical  meaning  of  the 
bending  stress  intensity  factors  and  K2  under  Sih's  defini- 
tion has  not,  to  the  author's  knowledge,  been  clarified  yet. 
First,  we  briefly  describe  the  process  of  Sih's  definition 
and  redefine  them  to  make  their  physical  meaning  clear. 

The  sum  of  the  normal  stresses  is, 
for  stretching: 

ax  + ay  = kj^  (p)  ^cosl  - k2(p)issin|-  (4.69) 

and  for  bending 
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°r  + c9  = 


12 (Mx  + M ) z 


= - — sin4] 

(3  + v ) h /r  2 /r  2 


(4.70) 


JL 


where  and  k2  are  stress  intensity  factors  for  plate  exten- 
sion and  and  K2  are  bending  stress  intensity  factors  defined 
by  Sih  such  that  both  bending  and  stretching  stress  intensity 
factors  have  the  identical  character.  These  stress  intensity 
factors  are  represented  as  the  real  and  negative  imaginary 
parts  of  a complex  constant  as: 


k = k1  - ik2 


(4.71) 


K = Kl  - iK2 


(4.72) 


1.  0 

Introducing  the  complex  variable  z = re  , Eguations  (4.69) 
and  (4.70)  can  be  rewritten  as, 


for  stretching: 
ax  + ay  = Retk^)1*} 


(4.73) 


and  for  bending 


M + M = (1-+  V)h  Re{K(^)**} 
x y 3/2(3  + v)  2 


(4.74) 


The  relationship  between  the  stress  functions  and  the  left- 
hand  members  of  Equations  (4.73)  and  (4.74)  are 
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for  extension: 


a + a = 4Re{4>  ' (z)  } 

a y s 


(4.75) 


and  for  bending 


Mx  + My  = -4D(1  + v)Re{<tfa' (z)  } (4.76) 

where  $g(z)  and  $b(z)  are  stress  function  for  plate  stretching 

and  for  plate  bending,  respectively.  In  the  vicinity  of  a 

crack  tip,  the  dominant  behavior  of  $'(z)  takes  the  form  of 

1/2 

a complex  constant  divided  by  z . Therefore,  equating  the 
results  of  Equations  (4.73)  and  (4.75),  the  stress  intensity 
factors  can  be  written,  for 
stretching: 


k » 2/2  lim  zS'  (z) 
z-0  s 

and  for  bending 
„ _ 12/5D(3  + v).  . 


lim  z 4.  ' (z) 
z-0  b 


(4.77) 


(4.78) 


Equation  (4.78)  is  the  definition  of  the  bending  stress 
intensity  factors  proposed  by  Sih  et  al.  and  Kj  are 

called  the  bending  stress  intensity  factor  and  the  shearing 
stress  intensity  factor,  respectively. 

The  Equation  (4.78) involves  the  material  constant  v 
explicitly  and  so  the  above  definition  can  be  applied  only 
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for  the  plate  of  isotropic  material.  The  present  author 
defines  the  stress  intensity  factors  for  flexural  problems 
in  the  following  manner. 


K.  = lim  - — r1*  M (r,  0) 
rn-0  h y 

K,  = lim  ^~r3/z  V (r,  0) 
* r-0  ti  y 

9M 

where  Vy  = Qy  + -g* 


(4.79) 


The  stress  intensity  factors  thus  defined  can  be  shown 
to  coincide  with  those  defined  by  Sih  in  isotropic  problems. 
Their  physical  meaning  is  apparent  i.e.,  represents  the 
intensity  of  bending  moments,  whereas  K2  corresponds  to 
statically  equivalent  stress  resultants  in  the  Kirchhoff's 
sense.  This  definition  can  be  immediately  applied  to  the 
flexural  problem  of  the  plate  made  of  anisotropic  materials. 
From  the  definition  (4.79),  together  with  equations  (4.26) 
and  (4.27),  the  bending  and  vertical  force  stress  intensity 
factor  K.^  and  can  be  written  in  terms  of  B's  as 

„ _ 6/2D(3  + v)„ 

K1 75  »i 


. 6/2D(3  + v) 

2 h2  P2N+1 


(4.80) 


Since  the  field  parameters  B's  are  derived  from  the  general- 
ized nodal  displacement  vector  q using  Equation  (4.55),  and 
K2  are  expressed  in  the  form 


h Bq 


(4.81) 


For  the  crack  tip  superelement  shown  in  Figure  32 , the  matrix 
B is  given  in  Table  30  for  E=l,  v=0.3. 

4.11  Numerical  Examples 

Several  performance  tests  have  been  carried  out  in  order 
to  investigate  the  accuracy  of  the  stress  intensity  factor 
solution  and  the  effect  of  the  size  of  the  crack  tip  super- 
element.  In  the  first  example,  the  crack  tip  superelement  has 
been  applied  to  the  problems  for  which  exact  solutions  exist. 
As  the  second  example  an  analysis  has  been  made  of  the 
bending  of  the  center  cracked  plate  of  finite  dimension,  for 
which  the  finite  element  solution  was  obtained  by  Wilson  and 
Thompson  [109]  by  using  the  conventional  elements. 

Example  1 

The  plate  geometry  is  a certain  cracked  plate  of  infinite 
dimensions.  The  through- the- thickness  crack  of  a total  length 
of  2a  lies  along  the  x-axis.  Utilizing  the  symmetries  in  the 
problem,  one  half  of  the  plate  is  concerned.  A nine-node 
square  shape  superelement  occupies  a crack  tip  region  as  shown 
in  Figure  33.  The  total  length  of  the  crack  is  varied  from 
2e  to  20e  where  e is  the  length  of  the  crack  within  the  super- 
element  and  the  side  of  the  element  is  equal  to  2e. 


* 
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Three  different  loading  conditions,  i.e.  purely  cylindri- 


cal  bending,  pure  bending  and  twisting  are  considered.  The 
closed  form  analytical  solution  for  the  infinite  plate  with 
an  elliptical  hole  is  readily  available  [113]  in  the  form  of 
a stress  function  under  above  loading  conditions.  The  through- 
the-thickness  crack  problem  is  a limiting  case  of  this  problem. 
Using  these  analytical  solutions  of  closed  form,  the  deflec- 
tions w,  w , w are  calculated  at  all  points  on  which  the 

> x i y 

nodal  points  of  the  superelement  are  located.  With  these 
corresponding  values  as  prescribed  displacements,  the  stress 
intensity  factor  solutions  are  computed  from  Equation  (4.81) 
for  several  ratios  a/e  and  compared  with  the  exact  stress 
intensity  factor  solutions  given  in  the  handbook  [33] . 

Case  1 Pure  Cylindrical  Bending  of  an  Infinite  Plate  Contain- 
ing a Finite  Crack 

This  problem  is  illustrated  in  Figure  33(a).  The  uni- 
form bending  moment  Mq  per  unit  length  is  applied  at  infinity 
along  the  edges  which  are  parallel  to  the  x-axis  on  which  the 
crack  is  located  >..*• 


M = M. 

y o 


The  exact  stress  intensity  factor  solutions  for  this  problem 
are 
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where  h is  the  plate  thickness. 

The  computational  results  are  tabulated  in  Table  25  for 
several  ratios  a/e.  It  is  observed  that  the  best  estimation 
for  the  value  is  obtained  when  the  superelement  contains 
the  entire  crack.  In  that  case  the  element  gives  an  error 
of  only  0.18  percent.  However,  when  the  superelement  is 
taken  as  small  as  one- tenth  of  the  size  of  the  length  of  the 
crack  the  solution  still  shows  an  excellent  accuracy. 

The  error  is  1.22%.  The  values  of  K2  solution  are  negligibly 
small  for  the  entire  yange  of  a/e. 


Case  2 Pure  Bending  of  an  Infinite  Plate  Containing  a 


Finite  Crack 


Figure  33(b)  illustates  this  example.  Uniform  bending 
moment  of  Mq  per  unit  length  is  applied  all  around  the 
boundary  of  the  plate  at  infinity 


The  exact  stress  intensity  factor  solutions  are  the  same  as 
those  in  the  previous  case. 

6M0  k 

K.  = — 7T-  ar  , K = 0 
1 ti  z 


The  results  are  independent  of  the  moment  M^.  The  finite 
element  solutions  are  shown  in  Table  25.  The  nature  of  the 
solutions  are  completely  the  same  as  that  in  the  case 
of  cylindrical  bending.  The  accuracy  of  the  K2  solutions 
appears  to  be  slightly  worse  but  it  is  still  of  the  same 
order. 

Case  3 Pure  Twisting  of  an  Infinite  Plate  Containing  a 
Finite  Crack 

Uniform  pure  twisting  moment  of  Mq  per  unit  length  is 
applied  all  around  the  contour  of  the  plate  at  infinity  as 
shown  in  Figure  33(c). 

M = M = 0 
x y 

Mxy  = M0 

The  exact  stress  intensity  factor  solutions  are 
= 0 

6M0  *5 

K2  = ~Ta 

z hi 

The  stress  intensity  factor  solutions  are  given  in  Table  25. 
It  is  seen  that  the  error  of  K2  solution  is  only  0.61%  when 
the  size  of  the  superelement  is  one-tenth  of  the  crack 
length  while  it  has  2.28  percent  error  when  the  superelement 


contains  the  entire  crack.  This  result  is  opposite  to  those 


in  the  previous  two  cases.  The  values  of  the  solution  is 
negligibly  small  in  every  a/e. 


It  can  be  seen  from  the  above  results  that  the  accuracy 
of  the  stress  intensity  factor  solution  is  excellent  and  it 
is  not  so  much  affected  by  the  size  of  that  crack  tip  super- 
element. The  size  of  the  superelement  may  be  of  the  same 
order  as  that  of  the  crack  length.  This  is  due  to  the 
accuracy  of  the  stress  and  displacement  field  in  the  interior 
of  the  superelement.  Further  improvement  of  the  solution  will 
be  achieved  by  increasing  both  the  number  of  the  parameter  3 ' s 
and  the  number  of  nodes. 

Example  2 

Consider  the  problem  of  pure  cylindrical  bending  of  a 
thin  plate  made  of  isotropic  materials  with  a finite  width 
and  centrally  located  through-the-thickness  crack,  whose 
geometry  is  shown  in  Figure  35.  Due  to  symmetry,  one- half 
of  the  plate  was  considered  in  the  finite  element  analysis. 

A number  of  cases  were  studied  with  different  ratios  of  2a/W, 
where  2a  was  the  total  crack  length  as  in  the  previous  example 
and  W,  the  width  of  the  plate.  In  each  case  only  a single 
nine-node  superelement  with  27  degrees  of  freedom  was  used 
[Figure  36].  Figure  37  presents  the  results  of  stress  inten- 
sity factor  versus  2a/W.  The  solid  line  in  this  figure 
represents  the  results  obtained  by  Wilson  and  Thompson  which 
were  computed  using  more  than  five-hundred  conventional  3- node 
triangular  plate  elements  based  on  assumed  displacement  finite 
element  model.  The  results  indicated  by  circles  are  those 

obtained  by  the  present  superelement. 

137 


It  is  seen  that  results  obtained  by  a single  crack  tip 
superelement  already  correlate  very  well  with  the  other 
results  which  were  estimated  to  be  accurate  to  within  2 
percent  according  to  their  paper.  In  this  simple  example 
problem,  it  is  shown  that  the  use  of  the  hybrid  superelement 
will  give  result  of  sufficient  accuracy  for  this  kind  of 
problem  even  when  only  a single  element  is  used. 

4.12  Conclusions  and  Recommendations 

A hybrid  finite  element  model  has  been  presented  for  the 
computation  of  bending  and  stress  resultant  stress  intensity 
factors  and  K2  respectively  in  the  flexural  problem  of  a 
plate  containing  a crack.  Prom  the  results  of  several  numeri- 
cal examples  the  following  conclusions  can  be  made: 

(1)  The  computation  time  required  for  the  generation 
of  a crack  element  stiffness  matrix  is  a fraction 
of  one  second  on  current  generation  large  computers. 
The  CPU  time  needed  for  the  entire  process  to  solve 
the  second  example  is  still  around  one  second. 

(2)  The  use  of  this  hybrid  crack  element  model  gives 
the  stress  intensity  factor  solutions  with  excel- 
lent accuracy. 

(3)  The  accuracy  of  the  stress  intensity  factor  solution 
is  not  considerably  affected  by  the  size  of  the 
crack  tip  superelement.  The  size  may  be  of  the 
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same  order  as  that  of  the  crack  length.  This 
indicates  that  the  very  crude  structure  models  can 
be  used  for  finite  element  analysis. 

(4)  It  is  to  be  noted  that  strickly  speaking  the  solu- 
tion to  the  flexural  problem  is  valid  only  when 
sufficient  extension  is  superimposed  such  that 
the  occurrence  of  the  overlapping  of  the  crack 
faces  on  the  compression  side  of  the  plate  can 
be  avoided. 

The  first  three  conclusions  are  valid  only  in  the  frame- 
work of  Poisson-Kirchhof f plate  theory.  A further  refinement 
of  this  superelement  is  possible  by  adopting  Reissner’s  sixth 
order  plate  theory.  In  his  theory,  the  deflection  w of  the 
plate  takes  the  form  of  a biharmonic  function  in  the  absence 
of  lateral  loads  and  can  be  represented  in  terms  of  two 
analytic  functions  of  complex  variable  z (z  = x+iy)  in 
exactly  the  same  form  as  that  in  the  classical  plate  theory. 

w = Re[z$(z)  + x(z)l  (4.82) 

The  average  rotations  8 and  8 about  y-  and  x-axes  respec- 

x y 

tively  are  given  by 


8 + i8  = -[4>(z)  + z$'(z)  - 2(n  - 1)  e <J>"  (z)  + x'(z)l 

x y 


+ i(n  - 1)  e2-^ 
3z 


(4.83) 
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where  n = - , e = h//TO 

and  h is  the  plate  thickness.  ip  is  the  stress  function  and 
satisfies  the  Helmholtz  equation. 

2 2 

- e V ip  = 0 


The  bending  moments  M , M , twisting  moment  M and  shear 

x y xy 

forces  Q and  Q can  be  expressed  in  terms  of  $ (z) , x(z)  and 
x y 

ip(z,z)  in  the  form 


Mx  + My  = DRe[*'(2)I 

Mx  - My  - 2iMxy  = tx"(z)  + z$"(z)  - 2(n-  1)  eV"  (z)  ] 

+ Sie2^ 

dz* 


Qx  - iQy  = -D[4$"  (z)  - 2if|] 


(4.84) 


The  functional  forms  of  the  complex  potential  4(z)  and  x(z) 
and  the  stress  function  ip(z,z)  are  determined  from  the  bound- 
ary conditions  of  the  problem.  The  procedure  to  obtain  the 
crack  element  stiffness  matrix  follows  completely  that  in  the 
classical  fourth  order  plate  theory. 
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SECTION  5 

BENDING  ANALYSIS  OF  BI-MATERIAL  PLATE 
WITH  THROUGH-THE-THICKNESS  CRACK  LOCATED 
ALONG  OR  NORMAL  TO  THE  INTERFACE,  ISOTROPIC 
PLATE  WITH  WEDGE-SHAPED  NOTCH  AND  ANISOTROPIC 
PLATE  WITH  THROUGH-THE-THICKNESS  CRACK 


1 


i 


5.1  Introduction 

In  this  section,  stress  and  displacement  fields  are  ob- 
tained for  bending  analyses  of  bi-material  plates  with  a crack 
located  along  or  normal  to  the  interface,  an  isotropic  plate 
with  a wedge-shaped  notch  and  an  anisotropic  plate  with  a 
through- the- thickness  crack.  The  procedure  to  construct  a 
superelement  for  above  cases  follows  completely  that  given 
in  Section  4. 

5.2  Bi-Material  Plate  with  a Through-the-Thickness  Crack 
Lying  along  the  Interface 

Figure  38  shows  a through- the- thickness  crack  lying  along 
the  interface  of  two  dissimilar  materials.  The  origin  of  the 
coordinates  is  placed  at  the  crack  tip.  The  positive  real 
axis  lies  along  the  line  of  bonding,  while  the  negative  real 
axis  lies  along  the  crack.  All  quantities  referred  to  the  re- 
gion y > 0 and  y < 0 are  designated  by  subscript  1 and  2,  re- 
spectively. 

On  the  basis  of  the  Poisson-Kirchhoff  theory  of  thin 
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plates  and  the  complex  variable  approach,  the  boundary  condi- 
tions can  be  stated  as  follows. 

(1)  Stress  free  condition  along  the  crack  surfaces  0 = ±ir 

(M  ) . = 0,  (V  ) . =0  at  0 » +7T 

Y 1 Y 1 (5.1) 

(M  )2  = 0,  (V  ) 2 = 0 at  6 = -7T 

(2)  Continuity  condition  of  the  bending  moment  M and  the 
Kirchhoff  shearing  force  V across  the  bonded  portion  of 
the  interface 

(V  1 = (V  2 at  9 - 0 

Y 1 y ^ (5.2) 

(V  1 = (Vy)  2 at  0 = 0 

(3)  Continuity  condition  of  the  deflection  w and  the  slope 

3w/3y  across  the  bonded  portion  of  the  interface 

(w)  .=  (w)  , at  0 = 0 

1 * (5.3) 

Ow/3y)  1=  (3w/3y)2  at  0 = 0 

All  field  variables  are  expressed  in  terms  of  a set  of 
analytic  functions  of  z,  i.e.,  <i>(z)  and  \{z)  . $(z)  and  x(z) 

in  both  materials  are  assumed  as 


< V 1 = <V  2 
(V  1 = <V  2 


(5.2) 


(w)  1=  (w)  2 


(5.3) 


...  X. 

$ (z)  = E aj1Jz  1 

1 x> 
$2(z)  = E a^z  i 


X.+l 

X1(z)  = E c ^ 1 z 

1 , ■ \ X.+l 

X2(z)  = E c(2x)z 


(5.4) 


where  a and  c2^  are  unknown  real  or  complex 

coefficients  and  can  be  evaluated  from  the  above  boundary  con- 
ditions. X^  are  real  or  complex  eigenparameters  which  are 
taken  to  be  the  same  in  both  regions. 

Inserting  Equation  (5.4)  into  the  boundary  conditions 

gives 
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) = o 


(X 

- 

. , iXir  . - -iXTT. 

y1)  (a1e  + a1e  ) 

+ ( X + 

1)  (c^e1^  + o1e"ii,T 

(X 

- 

. , -iXir  , - iXir. 

y 2)  (a2e  + a2e  ) 

+ ( X + 

1)  (c2e-iX''  + 

(X 

+ 

. , iXir  - -iXir. 

Ui>  (a^e  - axe  ) 

+ ( X + 

. . , iXiT  - -iXir 

1)  (c1e  - c±e 

(X 

+ 

y 2>  (a2e  ” a2e  ) 

+ ( X + 

. . , -iXir  - iXir 

1)  (c2e  - c2e 

y{  (x 

- y 1)  (a^  + a^)  + (X  + 

1) (cx  + 

• } 

= (X  - 

y2^  ^a2 

+ a2)  + (X  + 1)  (c2 

y{  (X 

+ yx)  (a1  - i^)  + (X  + 

1) (CL  - 

} 

= (X  + 

y2}  (a2 

- i2)  + (X  + 1)  ( c 2 

ai 

+ 

a.^  + + c.^  = a2  + a 

+ c _ + 
2 2 

52 

(X 

- 

1)  (ax  - a1  - a2  + a2) 

+ ( X + 

1)  (Cj^  “ - c2  + c 

) = o 
) = o 
) = 0 


(5.5) 


where  the  superscripts  are  omitted  and  y,  y^  and  y2  are 
defined  as 


D.d  - v.) 

Y = d2(1-"2T 


3 + vi 
yi  1 - v. 


(i  = 1,  2) 


(5.6) 


Equation  (5.5)  consists  of  eight  equations  with  eight  unknown 


constants  a^,  a^,  c^,  c^  (i  = 1,  2) 


By  eliminating  c^,  c^. 


c2  and  c2  from  Equation  (5.5),  one  obtains  the  following 
equation. 


1 + yLe 


i2Xir 


-1-V2e 


-i2Xir 


,,  i2Xir.  ...  -i2Xir, 
yy-L  (1  - e ) -y2  (1  - e ) 


0 

0 


0 

0 


0 

0 


0 

0 


1 + M1e 


-i2X7r 


-l-y2e 


i2X7r 


-i2Xir.  ..  i2Xir. 
yy-L  (1  - e ) -y2  (1  - e ) 
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The  above  equation  has  nontrivial  solutions  if  and  only  if 


the  coefficient  determinant  vanishes,  i.e. 


X = n 


-i2XTT 


(n  — 1,  2,  . ... ) 


V>i  (u 2 + Y) 

P2  + 1) 


(=  p) 


(5.8) 


(5.9) 


It  is  observed  that  the  eigenvalues  corresponding  to 
Equation  (5.9)  are  complex  and  are  expressed  in  the  form  of 


X = (n  - ±)  ± i< 
n 2 

where 


(n  = 1,  2, 


(5.10) 


< = - (1/2it)  logtv^  (u2  + y)  /y 2 (yp1  + 1)  } 


(5.11) 


The  deflection  and  stress  couples  are  expressed  as 


follows. 


For  X = 1,  2, 


in  material  1 
N Xn+1 


N x +1  X - y. 

w = E b r n [cos  (X  - 1)0 - ^cosU  + 1)9] 

n=l  n n X + 1 n 

N X -1  n 

M = -D.  E b X r n [2(1  + v.)cos(X„  - 1)9 

x 1 . n n 1 n 

n— l 

+ (1  - v. ) { (X  - 1)  cos  (X  -3)9  - (X  -y.)  cos  (X  - 1)6}] 
In  n n 1 n 

N X -1 

M = -D.  E b X r n [2(1  + v.)cosQ  - 1)9 
y 1 . n n in 


- (1  - v,){ (X_  - l)cos(X  - 3)9  - (X_  - y,)cos(X_  - 1)9}] 


Q = -4D.  E b X (X  - l)r  cos(X  -2)9 

x 1 , n n n n 

n=i 

N X -2 

Q = 4DX  E bnXn(Xn-l)r  n sin(Xn-2)9 
* n=l 

and  in  material  2 
N X +1 

w = a E b r n [cos  (X  -1)9  - {(X  - y ~)  / (X  +1)  }cos  (X  +1)9) 
, n n n 2 n n 

n=l 

N X -1 

M = -aD-  E b X r [2 (1  + v,) cos  (X  -1)9 
x 2 , n n 2 n 

n=l 

+ (1  - v,)  { U - 1)  cos  (X  - 3)9  - (X  - y,_)  cos  (X  - 1)9}] 

2 n n n 2 n 

N X -1 

M = -aD,  E bXr  n [ 2 (1  + v_) cos (X„  - 1) 9 

y 2 , n n 2 n 

n=l 

- (1  - v_)  { (X  - 1)  cos  (X  -3)9  - (X  - u. ) cos  (X  - 1)9}] 

2 n n n 1 n 

Mxy  = «D2a~V2)  E1Vnr  " [ “ D sin (Xn  - 3)  9 

n=l 

- ( X - y,)  sin  (X  -1)9] 
n z n 

N X -2 

Q = -4aD_  E b X (X  - 1) r cos (X  - 2)9 
x 2 t n n n n 

n— 1 

N X -2 

Q = 4aD_  E b X (X  - 1)  r n sin(Xn-2)9 
y 2 i n n n n 

n=i 

where 


a = (1  + y1)/  (1  + y2) 

For  X = (n-^)  ±ix  (n  = 1,  2,....) 

in  material  1 

N X +1  i (X  -1)9 
w = Re  E b r [e  n 

n 


i(Xn  + 1)0 


- { (X  - t)  / ( X + 1)  }e 


* 


N 


X -1  

M = -D.  Re  E bXr  n [2(l  + v.)e  n 

x 1 , n n 1 

n— i 


i (X  -1)0 


+ (1- v.){(X  -l)e 
1 n 


N 


i (X  -3)9 


- (Xn-t)e 
i (X  -1)6 


i (X  -1)0 


}] 


X -1 

M = -D,R  E bXr  n [2(l  + v.)e 

y i e„,  n n 1 

n— I 


i (X  -3)0  i (X  -1)0 

- (1  - v.)  { (Xn  - l)e  n - (X  - t)e  n }] 
in  n 

N X -1  i{X  -3)6 

“xy  = Dl(1-'’lII,E,V»r  t(Xn-l)e 

n=i 

i(X  -1)9 
- (Xn-t)e  n ] 

N X -2  i (X  -2)0 

Q = -4D, Re  E b X Q - 1)  r e 
x 1 n=1  n n n 


Qy  = 4D1ImJibnXn(Xn-l)r 


X_-2  i(X_-2)0 
e 


and  in  material  2 


N X +1  i ( X -1)  0 
w = -ytoRe  E pb  r n [e 
n=l 

N X -1  _,~ 

M = yojD.Re  E pb  X r [2(l  + V-)e 
x 1 , n n ^ 

n=l 


i ( X +1)0 

- ((xn- y2p)/(xn+ x)  }e  n 1 


+ (1  - v.)  { (Xn  - l)e 
2 n 


i (X  -3)0 


- (Xn  - y2P)  e 


N 


X -1  ... 

M = ya)D_Re  E pb  X r n [2(l  + v„)e 
y 2 n=1  n n z 

i(X  -3)0  ... 

- (1  - v2)  { (X  -De  n -(X  -y2p)e 


N 


i (X_-l) 9 

i ( X _— 1 ) 9 
i (X_-l) 0 

i(Xn-l)0 
i (X„  - 3)9 


}] 


}] 


n 


Mxv  = “YuD2(l  - v2)lm  E pbnXn[ (Xn  - l)e 
* n=l 

i (X  -1)9 

" (Xn"y2p)e  n 1 

N X -2  i(X  -2)0 

Q = 4yujD_Re  E pb  X (X  -l)r  e 

x 2 , n n n 

n— l 
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5.3  Bi-Material  Plate  with  a Through-the-Thickness 
Crack  Lying  Normal  to  the  Interface 
Figure  39  shows  a through-the-thickness  crack  lying 
normal  to  the  interface  of  two  dissimilar  materials.  The 
origin  of  the  coordinates  Is  placed  at  the  crack  tip.  The 
negative  real  axis  lies  along  the  crack,  while  the  imaginary 
axis  lies  along  the  line  of  bonding.  All  quantities  referred 
to  the  region  x < 0 and  x > 0 are  designated  by  subscript  1 and 
2,  respectively. 

On  the  basis  of  the  Poisson-Kirchhof f theory  of  thin 
plates  and  the  complex  variable  approach,  the  boundary  condi- 
tions can  be  stated  as  follows. 

(1)  Stress  free  condition  along  the  crack  surfaces  0 = ±tt 


(M  ) . =0 

y 1 

(V  ) . =0 

y 1 


at  0 = ±tt 


(5.13) 


(2)  Continuity  condition  of  the  bending  moment  and  the 

Kirchhoff  shearing  force  across  the  bonded  portion  of 
the  interface 


(Mx>  1 = <Mx>  2 
(Vl  " <Vx}2 


at  0 - 
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(5.14) 


LjtfT*  - 


(3)  Continuity  condition  of  the  deflection  w and  the  slope 


3w/3x  across  the  bonded  portion  of  the  interface 


(w)  1 = (w)  2 
(3w/3x)^  = (3w/3x) 2 


at  e = ±-2 


(5.15) 


(5-16) 


All  field  variables  are  expressed  in  terms  of  a set  of 

analytic  functions  of  z,  i.e.,  3>(z)  and  x(z).  ( z ) and  x(z) 

in  both  materials  are  assumed  as 

N , . . X . N , • , X . +1 

$x(z)  = Z a^1  z 1 X2(z)  = 2 ;z  1 

N ...  X.  N m X . +1  (5-16) 

<J>2(z)  = Z a2  z 1 X2(z)  = 2 c2  z 

where  a|^  , c|^  , a2^  and  c^1^  are  unknown  real  or  complex 
coefficients  and  can  be  evaluated  from  the  above  boundary 
conditions.  X^  are  real  or  complex  eigenparameters  which  are 
taken  to  be  the  same  in  both  regions. 

Symmetric  Loading  with  Respect  to  the  x-Axis 

In  this  case  a^)  and  c^^  are  real  constants  if  X^  are 
considered  to  be  real.  Substituting  Equation  (5.16)  into 
the  boundary  conditions  and  omitting  the  superscripts  yield 
the  following  equations. 

(e  +y1e  )a1~2Xa1~(e  - X)  a2  - (X  + 1)  c2  = 0 

2XaL  - (e"lXlT  + y1e~2lX7T)a1  - (X  - e"lX7r)a2  + (X  + l)c2  = 0 

(5.17) 

[ (Y y-L  + l)eiX7r  - (y  - l)u1e2iXlT]a1  + 2(y  - l)XaL  - (y2  + l)elXTTa2  = 0 


-iX-p- 


-2iXiT.  - 


2(y-l)Xa^  + [(yy^  + l)e  - (y-  Dy^  ]a1“(y2  + l)e 


-iXir 
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The  above  equations  have  nontrivial  solutions  if  and 
only  if  the  coefficient  determinant  vanishes,  i.e., 
sinAfr  = 0 

or  (5.18) 

41 2 (y  - 1)  (y  + y2)  + 2u1  (Y  + U2)  (Yy]_  + D cosXtt 

- y2 (y  - l)  (y  + u2)  + (yi^  + l)  (m2-yp1)  = 0 

All  unknown  constants  are  expressed  in  terms  of  a2(  = b) 
in  the  form  of 
a^  = f (X ) b 

= g ( X ) b (5.19) 

c 2 = h (X)  b 

where  b is  real  if  > is  considered  as  real  and 
f (X)  = (u2  + 1)  [ (YUX  + 1)  - (Y  - 1)  (Ux  + 2X)e“lX7r]/S(X) 
g (X)  = -(u2  + l)  [(Yl^  + 1)  U -y1e"2lX7r) 

+ (y  - 1)  (2X  - uL)  (X  - y^e"1^)  ] / (X  + l)S(X) 
h(X)  = (X  - e"lXlT)/(X  + 1)  - (y2  + 1)  (2XY(y1  + 1)  + u1(Y  - D (5.20) 
- [(Y  - 1)  (4X2  - y2)  + (Yyx  + 1)  ]e_lX7T 
+ y1(Yy1  + l)e"2iX7r}/(X  + l)S(X) 
s (X ) = (Yyx  + 1)  2 - 2yx  (Y  - 1)  (Yyx  + DcosXtt  + (Y  - 1)  2(yJ  - 4X2) 

The  resulting  displacement  and  stress  couples  are  in 
material  1: 

N X +1 

w = Re  E b r n [f_cos(X  - 1) 0 - f _sin (X  -1)0 
. n R n In 

n=l 

+ g_cos(X  + 1)  0 - g_sin (X  +1)0] 

R n In 
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{ [2(1  + vL)fR 


M = -D, Re  E X b r 

x 1 , n n 

n=l 

+ (1  - v1)  (Xn  + 1)  gR]  cos  (Xn  -1)0-  [2(l  + v1)f]. 

+ (1  - x>±)  (X^  + Dg-jJ  sin(Xn  - 1)  0 

+ (1  - v.)  U - 1)  [f 0cos (X  - 3)0  - fTsin(X  - 3)0] } 
In  R n i n 

N X -1 

M = -D.Re  E X b r n {[2(l  + v.)fn 
y 1 i n n l R 

n=l 

- (1  - v1)  (Xn  + l)gR]cos(Xn-  1)0  - [2(1  + vx)  fj 

- (1- v,)  (Xn  + l)g_]sin(X  - 1)0 

1 n I n 

- ( 1 - v, ) (X  - 1)  [f  cos (X  - 3)0  - fTsin(X  - 3)0]  } 

In  R n xn 

N X -1 

M = D,  (1  - v,  )Im  E b X r {{A  -1)  [f  sin(X  -3)0 
xy  x ^ ^ ii  n n rv  ri 

+ f _cos  ( X -3)0]  + (X  + 1)  [g  sin  (X  -1)0+  g cos  (X 
xn  n k n jl 

«x  - -4DlRe  ”1bnXn(Xn-1)'£Rc:os(Xn-2)e-£Isin<in-2 
n=l 

N 

Q = 4DnRe  E b X (X  -l)[f_sin(X  - 2)  0 + f cos  (X  - 2) 
1 , n n n R n I n 

1 n=l 

and  in  material  2: 

N X +1 

w = Re  E b r [cos  (X  - 1) 0 + hDcos (X  + 1) 0 ] 

, n n k n 

n=l 

N X — 1 

M = -D-Re  E X b r {2 (1  + v0) cos  (X„  - 1) 0 

x 2 n n 2 n 

n— 1 

+ (1  - v.)  [ (X  - l)cos(X  - 3)  0 + h_(X  + l)cos(X  - 1) 


L5 


N X -1 

M = -DnRe  E X b r { 2 (1 + v.) cos (X  -1)0 

y 2 n=l  n n 2 n 

-(l-v,)[(X  - 1)  cos  (X  - 3)  0 + h_U  + l)cos(\  - 1)  9]  } 

2 n n k n n 

N X -1 

M = D,(l  - v,)Re  E X br  [ U - 1)  sin  U - 3)  0 
xy  2 2 n_^  n n n n 

+ hR(Xn  + 1)  sin(Xn  - 1)  0] 

N X -2 

Q = -4D_Re  E b X (X  -l)r  cos  (X  -2)6 

x 2 , n n n n 

n— i 

N X -2 

Q = 4D0Re  E b X (X  -l)r  sin(X  -2)0  (5.21) 

y ' 2 n=l  n n n n 

where 

f (X)  = (u2  + 1)  [ (yy1  + 1)  - (y  - 1)  (y x + 2X)cosXir]/S  (X) 

fj  (X)  = (y2  + 1)  (y  - 1)  (y1  + 2 X ) sinXir/S  (X) 

gR(X)  = - (u2  + 1)  t (YP1  + 1)  (X  - u1cos2Xtt) 

+ (y-l)(2X-u1)(X  - u^cosXtt)  ] / ( X + 1)  S ( X ) 

gI(X)  = -v1  (u2  + 1)  [ (yyj^  + 1)  sin2Xi: 

+ (y  - 1)  (2X  - y ) sinXTT]/(X  + 1)S  (X) 

hR  (X)  = (X  - cosXtt)  / (X  + 1)  - (y2  + 1)  { 2Xy  (u1  + 1)  + y1(y-l) 

- [ (y  - 1)  ( 4X 2 - y2)  + (yyx  + 1)  IcosXtt 

+ y^^  (yy^  + 1)  cos2Xtt}/ (X  + 1)  S (X) 

S ( X ) = (yy^^  + 1)  2 - 2y1  (y  - 1)  (yy  ^ + 1)  cosXtt  + (y  - 1)  2 (y  2 - 4X2) 
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Antisymmetric  Loading 

The  constant  a ^ and  c^  are  pure  imaginary  if  X is  con- 
sidered as  real.  Imposing  boundary  conditions  gives  the  fol- 
lowing characteristic  equations. 


sinXir  = 0 


and  (5.22) 

4X  2 ( y - 1)  (y  2 + y)  +2y1(Y  + y2)  (yv1  + 1)  cosXtt  -y2(Y-l)(Y  + y2) 

+ (y\i1  + 1)  (y 2 - YU-^  = 0 

These  equations  are  identical  with  those  for  the  symmetric 
loading. 

All  unknown  constants  are  expressed  in  terms  of  a2(=ib) 
in  the  form  of 
a.^  = if  (X)  b 

C;L  = ig  (X)  b (5.23) 

c2  = ih(X)b 

where  b is  real  if  X is  considered  as  real  and 
f (X)  = (y2  + 1)  ((Y^!  + 1)  + (Y  - 1)  (2X  - y1)e'lX7T]/S(X) 
g (X)  = - (y 2 +D  (MYU-L  + 1)  + (Y  - 1)  (X  + yx)  (2X  - y1)e"iXlT 
+ y1(YU1  + l)e_2lX7T]/(X  + l)S(X) 

h(X)  = (X  + e‘lXTr)/(X  + 1)  - (y2  + 1)  {2XY(yx  + 1)  - y1(Y  - 1)  (5.24) 
+ ((Y-l)  (4X2-y2)  + (Yy1  + l)]e"iXlT 
+ y1(YU1  + l)e_2iXlT}/(X  + l)S(X) 

S(X)  = (YU x + 1)  2 - 2yx(Y  - 1)  (Yyx  + DcosXtt  + (Y  - 1)  2(y2  - 4X2) 

The  resulting  displacement  and  stress  couples  are  as 


follows. 


In  material  1: 


N X +1 

w = -Re  E b r n [f_sin(X  - 1) 0 + f_cos (X  - 1) 0 + g_sin(X  + 1) 0 
n R n In  R n 

n=l 

+ gj-cos  (Xn  + 1)  0 ] 

N X -1 

M = D.Re  E X b r n {2(1  + v.. ) [f-.sinQ  - 1)  0 + f Tcos U - 1)  0] 
x 1 , n n 1 R n In 


+ (1  - v. ) (X  - 1)  [f  sin(X  -3)0  + f cos (X  -3)0] 
in  R n in 


+ (1  - v. ) (X  + 1)  [g  sin ( X - 1)  0 + g cos(X  - 1)  0]  } 
x n x\  n x n 

N X -1 

= D.Re  E Xbr  {2(1  + v,)  [fDsin  ( X - 1) 0 + f Tcos ( X - 1) 0] 
y 1 i n n 1 R n In 


- (1  - v. ) (X  - 1)  [f  sin(X  -3)0  + f cos(X  -3)0] 
in  R n in 


- (1  - vx)  (Xn  + 1)  [gRsin(Xn  -1)0  + gIcos(Xn  - 1)  0 ] } 

N X -1 

M = D.  (1-  v.)Re  E Xbr  n { (X  - 1)  [fpcos(Xn  -3)0 
xy  1 1 n=l  n n n R n 

- f].sin(Xn  -3)0]  + (Xn  + 1)  [gRcos  (Xn  - 1)  0 - g;[sin(Xn  -1)0]} 
N 

Qv  = 4D, Re  E b X (X  - 1)  [fDsin(Xn  - 2) 0 + fTcos(Xn  - 2) 0] 

x 1 - n n n R n I n 

n=l 

N 

Qy  = 4D^Re  E bn^n(^n  “ 1)  [fRc°s  (Xn  - 2)  0 - f].sin(Xn  - 2)  0] 


and  in  material  2 
N X +1 

w = -Re  E b r [sin (X„  - 1) 0 + h_sin (X„  + 1) 0 ] 

In  n i\  n 

i i 

N X -1 

M = D9Re  E b X r {2  (1  + v0)  sin(X  -1)0 
x z «.  n n z n 


+ (1  - v2)  [(Xn  - l)sin(Xn-  3)0  + (Xn  + l)hRsin(Xn-  1)  0]  } 


N X -1 

M = D_Re  E b X r {2 (1 + v,) sin ( X„  - 1) 9 
y 2 n=1  n n 2 n 

- (1  - v-)  [Q  - UsinQ  - 3)  0 + (Xn  + l)hpsin(X„  - 1)  0]  } 
z n n n k n 

N X -1 

Mxy  = °2(1  “ v2)Re  Z bnXnr  [ (Xn  " 1)cos(xn  “ 3>  9 

+ (Xn  + l)hRcos(Xn  - 1) 0] 

N X -2 

<3  = 4D,Re  Z b X ( X - 1)  r n sin(X  -2)9 
x 2 n n n n 

n=l 

N X -2 

Q = 4D0Re  E b X (X  -l)r  n cos(X  -2)9  (5.25) 

y 2 i n n n n 

1 n=l 

where 

fR(X)  = (W2  + L)  t(Yyl  + 1}  + (Y  " 1}  (2X  " yl)cosXlTl/S(X) 
fI  ( X)  - -(u2  + l)  (Y-l)  (2X  - y^sinXir/SU) 

gR(X)  = - (u2  + 1)  [X  (yu1  + 1)  + (y  - 1)  (X  + y1)  (2X  - y^cosXir 
+ y^  (Yy^  + 1)  cos2Xir  ]/ (X  + 1)  S ( X) 
gj(X)  = (y2  + 1)  [ (y  - D (X  + yx)  (2X  - y ^ ) sinXir 
+ y.  (Yy,  + 1)  sin2Xir]/(X  + 1)S  (X) 
hR  (X)  = (X  + cosXtt)  / (X  + 1)  - (y2  + l){2XY(y1  + l)  - y^  (y  " D 
+ [ (y  - 1)  (4X2  - y2)  + (Yyx  + D ]cosXtt 
+ y^  (Yy^.  + D cos2Xtt}/ (X  + 1)  S (X) 
s (X)  = (Yyi  + 1)2-  2y1(Y-l)  (yy1  + l)cosXTT+  (Y-l)2(yJ-4X2) 


*»*■  ****■-  - -*7J  — 
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Consider  the  wedge-shaped  notch  shown  in  Figure  40.  The 
origin  of  the  coordinates  is  located  at  the  tip  of  the  notch 
and  the  notch  surfaces  are  situated  at  0 = ±ot  in  polar  coor- 
dinates . 

On  the  basis  of  the  Poisson-Kirchhof f . theory  of  thin 
plates  and  the  complex  variable  approach,  all  field  variables 


w,  M , M , M , Q and  Q are  expressed  in  terms  of  a set  of 
x y xy  x y 


analytic  functions  4>(z)  and  x(z),  which  are  assumed  as 
N X 

$(z)  = Z a z n 
n=l  n 
N 


X (z)  = l c z 
n=l  n 


V1 


(5.26) 


Stress  free  conditions  along  the  notch  surfaces  0 = ±ct 
can  be  expressed  in  a single  complex  equation  as 

(z)  + z4> ' ( z)  + x*  (z)  =0  (5.27) 

where  < = - (3  + v)/(l-v) 
v is  the  Poisson's  ratio. 

Inserting  Equation  (5.26)  into  Equation  (5.27),  one  ob- 
tains the  following  simultaneous  equations. 

Kae2lXa  + aX  e2ia  + c(X  + l)  = 0 

-2iXa  -n  -2ia  - (5*28) 

<ae  2lAa+aXne  2 +c(X  + l)  = 0 

Eliminating  c from  the  above  equations  gives  the  follow- 
ing characteristic  equation, 

<asin2Xa  + aXsin2a  = 0 (5.29) 

For  a symmetric  loading,  a and  c are  real  if  X is  treat- 


ed as  real.  The  characteristic  equation  (5.29)  yields 


icsin2Aa  + Asin2a  = 0 


(5.30) 


Putting  a = b (b  is  real.)  and  substituting  a = b into  Equa- 
tion (5.28)  gives 

c = -b  (<cos2Aa  + Acos2a) / ( A + 1)  (5.31) 


The  resulting  displacements  and  stress  couples  are 
N A +1 

w = Re  E b r [cos  (A  -1)0 
n=l  n n 

- t — t(kcos2A  a + A cos2a)cos(A  +1)0] 

A + l n n n 

n 

N A -1 

M = -DRe  E b A r [2(l  + v)cos(A  -1)0 

x , n n n 

n=l 

+ (1—  v)  {(A  - 1)  cos  (A  -3)0  - (kcos2A  a 

n n n 

+ A cos2a)cos(A  -1)0}] 
n n 

N A -1 

M = -DRe  E b A r [2(l  + v)cos(A  -1)0 

y n=l  n n n 

- (1—  v)  {(A  - 1)  cos  (A  -3)0  - (kcos2A  a 
n n n 


+ A cos2a)cos(A  -1)0}] 
n n 

N A -1 

M = D(l-v)Re  E b A r [(A  -l)sin(A  -3)0 
xy  . n n n n 


- (kcos2A  a + A cos2a)sin(A  -1)0] 
n n n 

N A -2 

Q = -4DRe  E b A (A  -l)r  cos(A  -2)0 
x n=1  n n n n 

N A -2 

Q = 4 DRe  E b A (A  -l)r  sin(A  -2)0 
Y n n n n 


(5.32) 


For  an  antisymmetric  loading,  a and  c are  pure  imaginary 
if  A is  considered  as  real.  Putting  a = ib  (b  is  real.)  and 


4 


substituting  a = ib  into  Equation  (5.28)  gives 
c = b (<cos2Aa  - Acos2a)  / (A  + 1) 


(5.33) 


The  resulting  displacement  and  stress  fields  are 
N X +1 

w = -Re  Z b r [sin(A  -1)0 
n=l  n n 

l 

- t — tt(<cos2A  a-A  cos2a)sin(A  +1)6] 

An  + 1 n n n 

N A -1 

14  = DRe  Z b A r (2(l  + v)sin(A  -1)0 

x . n n n 

n=l 

+ (1  - v)  [ (A  - 1)  sin  (A  -3)0  + (kcos2A  a 
n n n 


- A cos2a)  sin  (A  -1)0]} 

n n 

N A -1 

M = DRe  Z b A r n (2(l  + v)sin(A  -1)0 

y n=l  n n n 

- (l-v)l(X  -l)sin(A  -3)0+(kcos2A  a 

n n n 


- A cos2a)  sin  (A  -1)0]} 
n n 

N A -1 

Mxy  = D ( 1 - v)  Re  Z bnAnr  n t(An~  l)cos(An-  3)0 


+ (kcos2A  a-A  cos2a)cos(A  -1)0] 
n n n 

N A -2 

Q = 4D  Z b A (A  - 1) r n sin(A  -2)0 
x n n n n 

n=l 

N A -2 

Qy  = 4D  Z bnAn(An-  l)r  n cos(An-2)0  (5.34) 
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ar: 


Consider  a thin  plate  made  of  an  anisotropic  material 
which  possesses  a plane  of  elastic  symmetry  parallel  to  the 
midplane  of  the  plate  [Figure  28] . In  this  case,  the  consti 
tutive  relations  can  be  written  in  the  following  matrix  form 
e = Ca  (5.3 

where 


S * {sx  £y 

£z  Y 

yz  Y 

Y 

zx  'xy 

} 

2 * {ox  °y 

az  T 

yz  T 

zx  Txy 

} 

"cll 

C12 

C13 

0 

0 

C16 

C12 

C22 

C23 

0 

0 

C26 

c = 

C13 

C23 

C33 

0 

0 

C36 

0 

0 

0 

C44 

C45 

0 

0 

0 

0 

C45 

C55 

0 

_ C16 

C26 

C36 

0 

0 

C66 

The  above  elastic  constants  c.^,  ci2'  c13'  c22'  c23'  C33'  c4- 


c . _ , c _ _ and  c,,  can  be  expressed  in  terms  of  the  usual  engi' 

4 D D 0 O O 


neering  constants  according  to  the  following  formulae: 


cn  = 1/E  , c..  = -v  /E  = -v  /E 

11  ' x 12  xy'  y yx'  x 


c, . = -v  /E  = -v  /E  , c„ „ = 1/E 
13  xz'  z zx/  x 22  y 


c__  = -V  /E  = -V  /E  , c,-  = 1/E 

23  yz  z zy  y 33  z 


(5.3* 


c . . = 1/G  , ccc  = 1/G  , c,,  = 1/G 

44  ' yz  ' 55  ' zx  66  xy 


where  Ex  is  the  Young's  modulus  along  the  x-axis,  etc.  G^ 


is  the  shear  modulus  in  the  x-y  plane,  etc.  v is  the 


Poisson's  ratio  indicating  the  contraction  along  the  x-axis 
during  the  expansion  along  the  y-axis,  etc. 

Inverting  the  compliance  matrix  C,  the  Hooke's  law  can 
be  rewritten  in  the  following  matrix  form: 

a = De  (5.37) 

where  D is  the  stiffness  matrix  written  in  the  form  of 


11 

d12 

d13 

0 

0 

d16 

12 

d22 

d23 

0 

0 

d26 

13 

d23 

d33 

0 

0 

d36 

0 

0 

0 

d44 

d45 

0 

0 

0 

0 

d45 

d55 

0 

16 

d26 

d36 

0 

0 

d66 

The  bending  moments,  twisting  moments  and  shear  forces 
due  to  the  transverse  deflection  w(x,  y)  can  be  expressed  as: 


Mx  _12(all 


^ 2 
9 w 


9x‘ 


+ a 


9 2w 


12 


+ 2a 


9y‘ 


16 


3 2 2 

m — h /,  9w  9w  . 

My  " I2(a12^1  + a22^T 


16 


-2 

_9_w_) 

9x9y 

2 

9 w } 
9x9y 


M = — 1-55  (a.  , 
xy  12  16 


- 2 
9 w 


+ a 


9x 


26 


ij£  + 2a  -£*-) 

2 + Za66,  ’ 


(5.38) 


9y 


>,3  - 3 „ 3 

« _ h 9 w , 9 w 

x ‘ 3ai6i^ 

.3 


9x9y 


+ (a  '+  2a66)-^,  + a 6^] 
12'  66  9x9y2  269y3 


1 W in  3 rv  3 

~ h.  d W , , , n * d W 

Qy  "I2ta16^„3  + a12  2a66 


3x~ 


. 2»  + 3a26 
9x  9y 


3 

9 w 


9x9y 


2 + a22 


93w 

3y3 


where 
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11 

(dlld33 

’ dL>/d33  ' 

al2  = 

(d12d33 

d13d23)/d33 

16 

{d16d33 

d13d36)/d33 

' a22 

= (d22d33 

d23)/d33 

26 

(d26d33 

d23d36)/d33 

' a66 

(d66d33 

- d?6>/d33 

The  fundamental  equation  for  the  bending  of  thin  aniso- 
tropic plates  are 

->4  -4  „4  ~4 

9 w , „ _ 9 w ,o/_  , . x 9 w ^ A _ 9w 

11«  4 16„  3„  12  66  » 2„  2 26-  .3 

9x  3x  9y  9x  9y  9x9y 


- 4 

. a - n 

a22  4 u 

^9y 


(5.39) 


The  general  solution  of  the  above  equation  (5.39)  depends 
on  the  roots  of  the  following  characteristic  equation. 


4 3 2 

a__s  + 4a_-s  + 2(a...  + 2a.c,)s  + 4a1cs  + an . = 0 (5.40) 

22  26  12  66  16  11 


On  the  basis  of  energy  considerations,  it  can  be  shown 
that  the  above  characteristic  equation  has  no  real  roots. 
These  roots  can  be  expressed  in  the  form: 

- a1  + i61 


s2  = a2  + i02 

s3  = «i  - iBj_  = s± 
S4  = a2  " iB2  = *2 


(5.41) 


where  a^,  a2,  8^  and  62  are  real  constants  and  it  is  always 
possible  to  arrange  that  6^  > 0 anu  82  > 0. 

The  general  solution  for  Equation  (5.39)  can  be  express- 
ed as: 
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(5.42) 


w(x,  y)  = 2Re[F1(z1)  + F2  ( z2>  ] 
where 

z^  = x + s^y  = x + ot^y  + iB-^y 
z2  = x + s2y  = x + a2y  + iB2Y 

Introducing  $(z1)  = dF^/dz^^  , 1'(z2)  = dF2/dz2  and  substi' 
tuting  the  displacement  function  w(x,  y)  into  Equation  (5.38) 
gives  the  general  expression  for  moments  and  shear  forces  in 
terms  of  the  two  analytic  functions  $(z^)  and  ^ ( z 2 ) in  the 
form: 

h3 

Mx  = -^-Re[Pl$*  (zx)  + q14',(22)l 
h3 

My  = — g-Re  [p2<t  ’ (z.^  + q2vi'  ’ ( z2)  ] 
h3 

Mxy  - ~Re[p3*'  (zx)  +q34',(z2)]  (5.43) 

h3 

Qx  = ~Re[Slp4<r  (zx)  + s2q4r*(z2)] 
h3 

Qy  = -g-P.e[p4*"  (zx)  + q44'"(z2)] 
where 


Pi 

= 

all 

+ 

2 

a12Sl 

+ 

2a16sl 

ql 

= 

all 

+ 

a12S2 

+ 

2a16S2 

P2 

= 

a12 

+ 

a 2 2 S 1 

+ 

2a26Sl 

q2 

= 

a12 

+ 

a22S2 

+ 

2a26s2 

P3 

= 

al6 

+ 

a26Sl 

+ 

2a66sl 

q3 

s 

a16 

+ 

a26S2 

+ 

2a66s2 

+ 2acc)s.  + a,..s* 


p.  a . . / s.  + 3a. c 


+ (a.. 


q4  = all/S2  + 3al6  + (a12  + 2a66)S2  + a26S2 

2 2 

Introducing  new  variables  and  ^ (z^  = z ^ 

F1^Z1^  and  F2^z2^  can  be  exPressed  as 


Fl(zl(cl))  j^-L [ ( e2  j+1  + l62N+2j+l)cl 


2 j+1 


+ (62j+2  + l62N+2j+2) *>1  1 


f2(22(C2))  * 2 MY^x!  + iY 


2 j+1 


2 j+1  ' 2N+2j+l'*2 

+ ( Y + iy  ) j-2 j + 2 

iy2j+2  iy2N+2j+2;  q2  1 


By  applying  the  Kirchhoff  stress  free  condition  along 
the  crack  surface,  all  y's  can  be  expressed  in  terms  of  B's 
Hence  all  the  field  variables  can  be  expressed  in  terms  of 
B's  only  as: 

w = 2 ? {B2i+1[ReU2j+1)  +AReU2j  + 1)  + DImU2j  + 1)] 
j=l  J 

+ 62j  + 2[Re(?2j  + 2)  +CRe(c2j  + 2)  + BIm ( ?2j  + 2)  ] 


- [lm(;2j+1)  - BRe(C2j  + 1)  + Clm(c2j  + 1)  ] 


2N+2 j+1 


- e0MJ.0^J.0[Im(c?j  + 2)  -DRe(<;2j+2)  + Aim  ( c2  j + 2)  ] } 


2N+2 j+2 


■ 2j  + 2 

'2  J 


3 N 

«x  = -§4  .£1<  C2j  + 1)  ( 2 j - l)B2j+1[Re(p1i;2j_3)  + ARe  (q1?  2 j"3) 
+ Dim  (qxi;2j-3)  ] + ( 2 j + 2)  ( 2 j ) B2  j + 2 [Re  (p^ 2 j_2) 


+ CRe(q.  £2j'2)  +BIm(q1c2j"2)  ] 


- (2j  + 1)  (2j  - l)B2N+2j+1[Im(Piqj  ')  -BRe(qiC2J  ') 
+ CIm(q1C2I,‘3)  1 - (2j  + 2)  ( 2 j ) 6 2N+2j  + 2 [Im  (p^* j"  2) 

- DReCqj^?^"2)  + Alm(q1?2:i"2)  ] } 

3 N 

= ~ .f1{(2j  + 1) (2j  - 1) e2j+ltRe(p2^1j~3)  + ARe(q2^2 


+ DIm(q2C2j_3)  I + ( 2 j + 2)  (2j)S2j  + 2[Re(p2qj  “) 

+ CRe(q2£2j"2)  + Blm(q2;2j"2) ] 

- ( 2 j + 1)  ( 2 j - 1)  B2N+2j+itIin(P2C2j~3)  ” BRe(q2C2j~3) 
+ Clm(q2i;2j_3)  ] - ( 2 j + 2)  (2j)  B2N+2 j + 2 [Im  (p2^2 j'2) 


2j-2. 


- DRe(q2C2j  2)  +AIm(q2C2j  2)  ] } 


= Z { ( 2 j + 1)  ( 2 j - 1)  S2i  + 1[Re(p3?2j"3)  +ARe(q3£ 

j=l  J 

+ Dlm(q3?2j"3) ] + (2j  + 2) (2j) B2j+2[Re(p3;2j"2) 

+ CRe(q3C2j-2)  + BIm  (q^  2 j'  2)  ] 

- ( 2 j + 1)  ( 2 j - 1)  B2N+2j+1[Im(p3?2j_3)  - BRe(q3;2j"3 
+ CIm(q3S2j“3) ] - (2j  + 2) ( 2 j ) 62N+2j+2[Im(p3C2j_2) 


- DRe(q3?2j"2)  + Alm(q3i;2j“2)  ] } 


25-5. 


= ~ E^{(2j  + 1)  ( 2 j - 1)  ( 2 j - 3)  62j  + 1[Re(SlP4^J  J) 


+ ARe(s0qj,?2j'5)  + DIm(s^q^:J  J)  ] 


2j“  5, 


+ ( 2 j + 2)  ( 2 j ) ( 2 j - 2)  S 


2j  + 2 [Re  ^S1P4^1  )+CRe(s2q4;2  ) 

+ Blm(s2q4;2j"4) ] 

- ( 2 j + 1)  (2j  - 1)  ( 2 j - 3)  S2N+2j+1  [Im(SlP4^:>-5, 

- BRe(s2q4£2j”5)  +CIm(s2q4^2j_5)  ] 

- ( 2 j + 2)  ( 2 j ) ( 2 j - 2)  62N+2j  + 2[Im(SlP4Ci:i”4)  ' (s^^3'4) 

+ AIm(s2q4?2^“4)  ] } 

h3  N o-|  c 

Qy  = fg  .IM  (2j  + 1)  ( 2 j - 1)  ( 2 j - 3)  e2j  + l[Re(p4Cl  } 

+ ARe(q4C2j"5)  +DIm(q4^j"5)  ] 

+ (2j  + 2)  (2j)  (2j  - 2)  32j  + 2 [Re(p4^j"4)  + CRe(q4C2j"4) 

+ Blm(q4^j"4)  ] - (2j  + 1)  (2j  - 1)  { 2 j - 3)  B2N+2j+1  [InKp^3"5) 

- BRe(q4S2j"5)  + CIm(q4C2j'5)  ] 

- ( 2 j + 2)  ( 2 j ) ( 2 j - 2)  S2N+2j  + 2 [Im(p4^3-4)  . DRe  (q^D"4) 

+ Alm(q4^j"4)  ] } (5.45) 

where 

A=  (f4g1  - f2g3)/(f 2g4  - f4g2) 

B = (fLf4  - f2f3)/(f2g4  - f4g2) 

C = (f3g2  - f xg4) / (f 2g4  - f 4g2) 

D = (g,g,  - g,g7)/(f-g,  - f .g.) 
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243  POLYNOMIAL  BASIS  STRESS  ASSUMPTION  FOR  A TYPE  B CRACK  ELEMENT 


24B  POLYNOMIAL  BASIS  STRESS  ASSUMPTION  FOR  A TYPE  A CRACK  ELEMENT 


456  POLYNOMIAL  BASIS  STRESS  ASSUMPTION  FOR  A TYPE  B CRACK  ELEMENT 


TABLE  6 (CONCLUDED) 


nxy 


[P]  (B 


51 B POLYNOMIAL  BASIS  STRESS  ASSUMPTION  FOR  A TYPE  A CRACK  ELEMENT 


TABLE  8 (CONCLUDED) 


-2Hzx+nz 


TABLE  9 (CONCLUDED) 


TABLE 


TABLE  12 


DISPLACEMENTS  RELATED  TO  SINGULAR  STRESS  TERMS 
Us  = (u,  V,  w}  = [Ls]{Kr  Kir  KIXI} 
where 


[Ls] 


|^{(S-8v)cos^- cos^}  (9-8v)  sin^  + sin^p}  0 


sir 


|{  (7-8v)sin§-  sin-y-} 


(-3+8v)  cos^  - cos-^-}  0 


TABLE  13 

LINEAR  ISOTROPIC  COMPLIANCE  MATRIX 


FOR 

THREE-DIMENSIONAL  CRACK 

ELEMENT 

e = {e  » e . 

£ f Y / y * 

Y } 

x y 

z 'yz  ' zx 

' xy 

= c{q  . a 

/ 0 $ T f T » 

T } 

- 

x y 

z yz  zx 

xy 

• 

where 

1/E 
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1/E  -v/E 
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-v/E 

-v/E  1/E 
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C = 
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0 0 

2 (1+v) /E 

0 

0 

9 

0 0 

0 

2 ( 1+v) /E 
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0 0 

0 

0 

2 (1+v) 

E;  ELASTIC  MODULUS 


v;  POISSON'S  RATIO 


TABLE  14 


EVALUATION  OF  DIFFERENT  TYPES  OF  ELEMENT  AND  SUPERELEMENT 
Type  A 


No.  of  Nodes 

No.  of  8's 

8 

24 

8 

30 

12 

30 

16 

45 

KI 

.94250 

.94241 

.98944 

1.0031 

KII 

1.0325 

1.0287 

1.0015 

.97760 

KIII 

1.0574 

1.0754 

.99874 

.99988 

Type  B 

No.  of  Nodes 

8 

12 

12 

16 

No.  of  B's 

24 

45 

51 

51 

KI 

1.0126 

.99841 

1.0040 

1.0040 

KII 

1.0279 

.99461 

1.0012 

1.0012 

KIII 

1.0019 

.99896 

1.0029 

. 1.0029 

Half  Element 

(Type  A + Type  B) 

Total  No.  of 

Nodes 

12 

20 

No.  of  Nodes 

8 (Type 

A)  +8 ( Type  B) 

12  (Type  A) +12 (Type  B! 

No.  of  8's 

24 

24 

30 

45 

KI 

l 

.98245 

.99608 

KII 

V 

1 

..0136 

.99647 

KIII 

] 

-.0257 

.99886 
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<1 


198 


0257  .99886 
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TABLE  20 

CYLINDRICAL  ROD  WITH  AXISYMMETRICALLY 
LOCATED  PENNY  SHAPED  CRACK 
UNDER  UNIFORM  TENSION 


An^rle  of  One  Sector 

9 = 10° 

0 

tr> 

II 

CD 

396  d.o.f. 

Kz  = 1.583 

Kt  = 1.556 

52  Element  Grid 

(1.20%  Error) 

(2.87%  Error) 

810  d.o.f. 

Kz  = 1.616 

Kt  = 1.589 

116  Element  Grid 

(0.868%  Error) 

(0.805%  Error) 

Benthem's  Solution 

Kj  = 1.602 

Benthem's  Solution 


1.602 


TABLE  21 

DISTRIBUTION  OF  STRESS  INTENSITY  FACTOR  ALONG 
CRACK  FRONT  FOR  EMBEDDED  PENNY-SHAPED  CRACK 
IN  CUBE  UNDER  UNIFORM  TENSION 


Angle  from 
the  X-axis 

9 (degree) 

KI 

KI/2a/a/7r 

m 

1.578 

0.9889 

15  - 30 

1.578 

0.9889 

30  ~ 45 

1.578 

0.9889 

45  - 60 

1.578 

0.9889 

60  ~ 75 

1.578 

0.9889 

75  ~ 90 

1.578 

0.9889 

TABLE  22 

DISTRIBUTION  OF  STRESS  INTENSITY  FACTOR  ALONG 
CRACK  FRONT  FOR  SEMI-CIRCULAR  SURFACE  CRACK 
IN  RECTANGULAR  PRISM  UNDER  UNIFORM  TENSION 


Angle  from 
the  X-axis 
9 (degree) 


- 15 


15  ~ 30 


30  ~ 45 


45  - 60 


60  ~ 75 


75  ~ 90 


KI/2a/a/7r 


1.820 


1.698 


1.652 


1.629 


1.618 


1.613 


1.141 


1.064 


1.035 


1.021 


1.013 


l.oii 





TABLE  23 


DISTRIBUTION  OF  STRESS  INTENSITY  FACTOR  ALONG 
CRACK  FRONT  FOR  QUARTER-CIRCULAR  CORNER  CRACK 
IN  RECTANGULAR  PRISM  UNDER  UNIFORM  TENSION 


TABLE  25 


STRESS  INTENSITY  FACTOR  SOLUTIONS  OF  AN  INFINITE  PLATE 
CONTAINING  A FINITE  CRACK  SUBJECT  TO  PURE  CYLINDRICAL 
BENDING,  PURE  BENDING  AND  TWISTING 


K /K 

computed7  exact 

e/a 

Pure  Cylindri- 
cal Bending 

Pure  Bending 

Pure  Twisting 

1/2 

0.9897 

0.9897 

0.9902 

1/3 

0.9887 

0.9887 

0.9921 

1/4 

0.9883 

0.9883 

0.9929 

1/5 

0.9881 

0.9881 

0.9932 

1/6 

0.9880 

0.9880 

0.9935 

1/7 

0.9880 

Q.9880 

0.9936 

1/8 

0.9879 

0.9879 

0.9937 

1/9 

0.9878 

0.9879 

0.9938 

1/10 

0.9878 

0.9878 

0.9939 

a;  Half  Length  of  the  Crack 

e;  Length  of  the  Crack  within  the  Superelement 
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TABLE  27  (CONTINUED) 
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TYPE  B ELEMENT 


FIG.  4 GENERAL  HEXAHEDRON  CRACK  ELEMENT 


6 


FIG.  5 GAPPING  OR  OVERLAPPING  BETWEEN  NEIGHBORING 

CRACK  ELEMENTS  WITH  A CURVED  CRACK  FRONT  EDGE 
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FIG.  b FAMILY  OF  THREE-DIMENSIONAL  HYBRID  STRESS  CRACK  ELEMENTS 


FIG.  9 TYPICAL  16-NODE  CRACK  FRONT  ELEMENT  AND  ITS 

RELATIVE  NODE  AND  NODAL  DISPLACEMENT  NUMBERING 
SCHEME 


FIG.  10  FRACTURE  TEST  SPECIMENS 


FIG.  11  FINITE  ELEMENT  BREAKDOWN  FOR  THE  ANALYSIS  OF  FRACTURE 
TEST  SPECIMENS 


NORMALIZED  STRESS  INTENSITY  FACTOR  Y 


35  ELEMENT  MESH  SUBDIVISION  WITH 
20-NODE  SINGULAR  HALF  ELEMENT 

35  ELEMENT  MESH  SUBDIVISION  WITH 
12- NODE  SINGULAR  HALF  ELEMENT 

85  ELEMENT  MESH  SUBDIVISION  WITH 
12-NODE  SINGULAR  HALF  ELEMENT 


E = 1 psi,  v = .3, 

b/a  = 2,  L/a  =1,  a = 2" , 2t  = 1" 
= Ycra*5 


MIDPLANE 


.25 

2 y/t 

DISTANCE  FROM  THE  MIDPLANE 

FIG,  12  SINGLE  EDGE  CRACK  SPECIMEN 


SURFACE 


NORMALIZED  STRESS  INTENSITY  FACTOR 


i TWO-DIMENSIONAL  SOLUTION 


35  ELEMENT  MESH  SUBDIVISION  WITH 
20-NODE  SINGULAR  HALF  ELEMENT 

35  ELEMENT  MESH  SUBDIVISION  WITH 
12-NODE  SINGULAR  HALF  ELEMENT 

85  ELEMENT  MESH  SUBDIVISION  WITH 
12-NODE  SINGULAR  HALF  ELEMENT 


E = 1 psi,  v = 


b/a  = 2 , L/a 
Kx  = Ycra*5 


= 2,  a = 2",  2 1 = 1 " 


MIDPLANE 

y/t 

DISTANCE  FROM  THE  MIDPLANE 

SURFACE 

FIG, 

13  CENTER  CRACK  S°ECIMEN 

NORMAL  I 


TWO-DIMENSIONAL  SOLUTION 


s 

A 


0 

MIDPLANE 


8 

A 


i a g 

AAA 


O 35  ELEMENT  MESH  SUBDIVISION  WITH 
20-NODE  SINGULAR  HALF  ELEMENT 

A 35  ELEMENT  MESH  SUBDIVISION  WITH 
12-NODE  SINGULAR  HALF  ELEMENT 

□ 85  ELEMENT  MESH  SUBDIVISION  WITH 

12-NODE  SINGULAR  HALF  ELEMENT 


E = 1 psi,  v = . 3 , 

b/a  = 2,  L/a  = 2,  a = 2" , 2t  = 1" 
Kx  = Yaa*5 


.2 

y/t  SURFACE 

DISTANCE  FROM  THE  MIDPLANE 


PIG.  14  DOUBLE  EDGE  CRACK  SPECIMEN 


H = 0.6W 
h = 0.2  7 5W 
D = 0.2  5 W 
c = 0.2  5W 
2 t = 0.5  0W 


FIG.  15  STANDARD  COMPACT  TENSION  SPECIMEN 

(astm  standard  e-399-72) 


FINITE  ELEMENT  BREAKDOWN  FOR  THE  ANALYSIS  OF  STANDARD 
COMPACT  TENSION  SPECIMENS 


NORMALIZED  STRESS  INTENSITY  FACTOR 


TWO-DIMENSIONAL  SOLUTION 


1.0 


□ 


YAMAMOTO'S  SOLUTION 


O 12-NODE  SINGULAR  HALF  ELEMENT 
COMPOSED  OF  TWO  8-NODE  BASIC 
SINGULAR  ELEMENTS 


05 


a/W  = 5/8,  v = .3, 


Y 


, T(2W  + a) 
I7  t (W  - a)  V2 


0 


5 


MIDPLANE 


y/t 


SURFACE 


DISTANCE  FROM  THE  MIDPLANE 


PIG.  18  STRESS  INTENSITY  FACTOR  FOR  COMPACT 
TENSION  SPECIMEN 


ITE  ELEMENT  BREAKDOWN  FOR  THE  ANALYSIS  OF  A PENNY-SHAPED  CRACK 
INDR  I CAL  ROD 


SINGULAR  ELEMENT 


K 

k 

FIG.  24  FINITE  ELEMENT  BREAKDOWN  FOR  THE  ANALYSIS  OF  A 
j PENNY-SHAPED  CRACK,,  A SEMI-CIRCULAR  SURFACE  CRACK 

AND  A QUARTER-CIRCULAR  CORNER  CRACK  IN  A 
1 RECTANGULAR  PRISM  IN  TENSION 


PRESENT  RESULT 
TRACEY  S RESULT  FOR  A 
PENNY-SHAPED  CRACK  IN  A 
CYLINDRICAL  ROD 
SNEDDON  S SOLUTION  FOR  A 
PENNY-SHAPED  CRACK  IN  AN 
INFINITE  SOLID 


ANGULAR  POSITION,  0-DEGREES 

FIG.  25  STRESS-INTENSITY  MAGNIFICATION  FACTOR  FOR  A 
PENNY-SHAPED  CRACK  IN  A CUBE  IN  TENSION 


STRESS-INTENSITY  MAGNIFICATION  FACTOR.,  KT/2o/aA 


SNEDDON 


O PRESENT  RESULT 

□ TRACEY'S  RESULT  FOR  A SEMI- 
CIRCULAR SURFACE  CRACK  IN  A ROD 
OF  SEMI-CIRCULAR  CROSS  SECTION 

A smith's  RESULT  FOR  A semi- 
circular SURFACE  CRACK  IN  A SEMI- 
INFINITE SOLID 


ANGULAR  POSITION.,  9-DEGREES 

26  STRESS-INTENSITY  MAGNIFICATION  FACTOR  FOR  A SEMI- 
CIRCULAR SURFACE  CRACK  IN  A RECTANGULAR  PRISM  IN 
TENSION 
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1 0.2 

CRACK  LENGTH/PLATE  WIDTH 

FIG.  37  BENDING  STRESS  INTENSITY  FACTOR  SOLUTIONS  FOR 
CENTER  CRACKED  PLATES  SUBJECTED  TO  CYLINDRICAL 
BENDING 
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